{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Secure Kaplan-Meier Survival Analysis Explained\n",
    "\n",
    "The MPyC demo [kmsurival.py](kmsurvival.py) implements privacy-preserving Kaplan-Meier [survival analysis](https://en.wikipedia.org/wiki/Survival_analysis), based on earlier work by Meilof Veeningen. The demo is built using the Python package [lifelines](https://pypi.org/project/lifelines/), which provides extensive support for survival\n",
    "analysis and includes several datasets. We use lifelines for plotting Kaplan-Meier survival curves, performing logrank tests to compare survival curves, and printing survival tables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# use pip (or, conda) to make sure lifelines package is installed:\n",
    "!pip -q install lifelines  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import os, functools\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import lifelines.statistics\n",
    "from mpyc.runtime import mpc\n",
    "mpc.logging(False)\n",
    "from kmsurvival import fit_plot, events_to_table, events_from_table, logrank_test, aggregate, agg_logrank_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actual use of the lifelines package is hidden mostly inside the [kmsurival.py](kmsurvival.py) demo, except for the function `lifelines.statistics.logrank_test()` used below to validate the results of the secure computations.\n",
    "\n",
    "## Kaplan-Meier Analysis\n",
    "\n",
    "We analyze the aml (\"acute myelogenous leukemia\") dataset, which is also used as a [small example on Wikipedia](https://en.wikipedia.org/wiki/Survival_analysis#Example:_Acute_myelogenous_leukemia_survival_data). The file `aml.csv` (copied from the R package `KMsurv`) contains the raw data for 23 patients. Status 1 stands for the event \"recurrence of aml cancer\" and status 0 means no event (\"censored\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style  type=\"text/css\" >\n",
       "</style><table id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783eb\" ><thead>    <tr>        <th class=\"col_heading level0 col0\" >observation</th>        <th class=\"col_heading level0 col1\" >time</th>        <th class=\"col_heading level0 col2\" >status</th>        <th class=\"col_heading level0 col3\" >group</th>    </tr></thead><tbody>\n",
       "                <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow0_col0\" class=\"data row0 col0\" >12</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow0_col1\" class=\"data row0 col1\" >5</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow0_col2\" class=\"data row0 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow0_col3\" class=\"data row0 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow1_col0\" class=\"data row1 col0\" >13</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow1_col1\" class=\"data row1 col1\" >5</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow1_col2\" class=\"data row1 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow1_col3\" class=\"data row1 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow2_col0\" class=\"data row2 col0\" >14</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow2_col1\" class=\"data row2 col1\" >8</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow2_col2\" class=\"data row2 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow2_col3\" class=\"data row2 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow3_col0\" class=\"data row3 col0\" >15</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow3_col1\" class=\"data row3 col1\" >8</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow3_col2\" class=\"data row3 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow3_col3\" class=\"data row3 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow4_col0\" class=\"data row4 col0\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow4_col1\" class=\"data row4 col1\" >9</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow4_col2\" class=\"data row4 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow4_col3\" class=\"data row4 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow5_col0\" class=\"data row5 col0\" >16</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow5_col1\" class=\"data row5 col1\" >12</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow5_col2\" class=\"data row5 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow5_col3\" class=\"data row5 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow6_col0\" class=\"data row6 col0\" >2</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow6_col1\" class=\"data row6 col1\" >13</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow6_col2\" class=\"data row6 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow6_col3\" class=\"data row6 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow7_col0\" class=\"data row7 col0\" >3</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow7_col1\" class=\"data row7 col1\" >13</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow7_col2\" class=\"data row7 col2\" >0</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow7_col3\" class=\"data row7 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow8_col0\" class=\"data row8 col0\" >17</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow8_col1\" class=\"data row8 col1\" >16</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow8_col2\" class=\"data row8 col2\" >0</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow8_col3\" class=\"data row8 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow9_col0\" class=\"data row9 col0\" >4</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow9_col1\" class=\"data row9 col1\" >18</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow9_col2\" class=\"data row9 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow9_col3\" class=\"data row9 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow10_col0\" class=\"data row10 col0\" >5</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow10_col1\" class=\"data row10 col1\" >23</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow10_col2\" class=\"data row10 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow10_col3\" class=\"data row10 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow11_col0\" class=\"data row11 col0\" >18</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow11_col1\" class=\"data row11 col1\" >23</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow11_col2\" class=\"data row11 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow11_col3\" class=\"data row11 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow12_col0\" class=\"data row12 col0\" >19</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow12_col1\" class=\"data row12 col1\" >27</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow12_col2\" class=\"data row12 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow12_col3\" class=\"data row12 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow13_col0\" class=\"data row13 col0\" >6</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow13_col1\" class=\"data row13 col1\" >28</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow13_col2\" class=\"data row13 col2\" >0</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow13_col3\" class=\"data row13 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow14_col0\" class=\"data row14 col0\" >20</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow14_col1\" class=\"data row14 col1\" >30</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow14_col2\" class=\"data row14 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow14_col3\" class=\"data row14 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow15_col0\" class=\"data row15 col0\" >7</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow15_col1\" class=\"data row15 col1\" >31</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow15_col2\" class=\"data row15 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow15_col3\" class=\"data row15 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow16_col0\" class=\"data row16 col0\" >21</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow16_col1\" class=\"data row16 col1\" >33</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow16_col2\" class=\"data row16 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow16_col3\" class=\"data row16 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow17_col0\" class=\"data row17 col0\" >8</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow17_col1\" class=\"data row17 col1\" >34</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow17_col2\" class=\"data row17 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow17_col3\" class=\"data row17 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow18_col0\" class=\"data row18 col0\" >22</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow18_col1\" class=\"data row18 col1\" >43</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow18_col2\" class=\"data row18 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow18_col3\" class=\"data row18 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow19_col0\" class=\"data row19 col0\" >9</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow19_col1\" class=\"data row19 col1\" >45</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow19_col2\" class=\"data row19 col2\" >0</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow19_col3\" class=\"data row19 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow20_col0\" class=\"data row20 col0\" >23</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow20_col1\" class=\"data row20 col1\" >45</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow20_col2\" class=\"data row20 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow20_col3\" class=\"data row20 col3\" >2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow21_col0\" class=\"data row21 col0\" >10</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow21_col1\" class=\"data row21 col1\" >48</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow21_col2\" class=\"data row21 col2\" >1</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow21_col3\" class=\"data row21 col3\" >1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                                <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow22_col0\" class=\"data row22 col0\" >11</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow22_col1\" class=\"data row22 col1\" >161</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow22_col2\" class=\"data row22 col2\" >0</td>\n",
       "                        <td id=\"T_ded87dde_d9fc_11e9_a697_e4b97af783ebrow22_col3\" class=\"data row22 col3\" >1</td>\n",
       "            </tr>\n",
       "    </tbody></table>"
      ],
      "text/plain": [
       "<pandas.io.formats.style.Styler at 0x1c352c0c848>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(os.path.join('data', 'surv', 'aml.csv')).rename(columns={'Unnamed: 0': 'observation', 'cens': 'status'})\n",
    "df.sort_values(['time', 'observation']).style.hide_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time is in weeks. The study compares the time until recurrence among two groups of patients. Patients in group 1 received maintenance chemotherapy, while patients in group 2 did not get any maintenance treatment. To plot the [Kaplan–Meier survival curve](https://en.wikipedia.org/wiki/Kaplan–Meier_estimator) for groups 1 and 2, we use the function `fit_plot()` as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEOCAYAAADc94MzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXwV9dn38c8FBMISFgUqsgh6iwIii2ERqIIIQrXaKq0ivZXiUu3jftfWCjcqlru2dtOn3rW4URX3rai4F4sPigqKCyAtspQgZZMlCGG9nj9mEg8nJ8kknOTMSb7v1+u8ODPzm5nrDMm5MjO/+V3m7oiIiMRNvUwHICIikooSlIiIxJISlIiIxJISlIiIxJISlIiIxJISlIiIxJISlNQ4M3Mz+49Mx1EWM7vRzO7NdByVZWYvmdmFadjOm2Z2cTpiEjkYSlCSkpmtNLNTE6bPM7PNZnZyJuMqj5l1DpPfB0nzW5vZbjNbGWU77v4/7p51X9DuPtrd/5LpOETSRQlKKhT+VX4XcLq7/z3T8UTQ1MyOS5g+H1hR3Tu1QLX8TplZg+rYbqbUts8j1UMJSsplZpcCvwVOc/e3E+Y/aWb/NrOtZjbHzHokLJtuZneb2WtmVmhmfzezI8rY/ulm9qGZbTOz1WZ2c8Ky4jOiC83sX2a20cwmRgj7ISDxUtcFwINJ+z3czJ42sw1mtsLMrkpYdrOZPZwwPdDM3jazLWb2kZkNTVj2pplNNbO5wA7gyBSf8WdmtiY8FkvNbHjCcfpFQruhZlaQML0yXPdj4Cszm2RmTyVt+w4zuzMhlovNrFEY63EJ7dqY2U4za2tmrczshfCzbw7fd4hwXDGz+uEl0M/Dz7PAzDom/F81SGhbcqnQzMab2Vwz+72ZfQncWl6M4fQZZrYwbPe2mR1f0TGV2kUJSspzOXArMNzd5yctewk4GmgLfADMSFo+Lly3NbAwxfJiXxEkkJbA6cDlZvadpDZDgGOA4cBkM+tWQdwPA+eFX6bdgDzg3eKF4VnO88BHQPtwu9eY2WnJGzKz9sCLwC+AQ4CfAE+bWZuEZv8JXBruZ1XS+scAVwD93D0POA1YWUH8icYSHJeWBIn3W2bWPNx2feD7wCOJK7j7LuCZcN1i3wf+7u7rCX7vHwCOADoBO4E/RoznunC73wKaAxMIEnMUA4DlBD8zU8qL0cz6AvcDPwIOBf4MzAyT78EeU8kSSlBSnhHAPOCT5AXufr+7F4ZfhjcDvcysRUKTF919Trh8InCimXVMsZ033f0Td9/v7h8DjwLJ97lucfed7v4RQVLpVUHcBcBS4FSCM6kHk5b3A9q4+xR33+3uy4F7gPNSbOsHwCx3nxXG+Bown+ALuth0d1/k7nvdfU/S+vuARkB3M8tx95Xu/nkF8Se6091Xh59/FcEfA8UJ/BRgh7vPS7HeIxz45X9+OA933+TuT7v7DncvBKZS+piX5WJgkrsv9cBH7r4p4rpfuPv/DY/TzvJiBC4B/uzu77r7vvDe2i5gIAd/TCVLKEFJeS4DugL3mpkVzwzPTG4LL/Ns4+u/XlsnrLu6+I27bwe+BA5P3oGZDTCz2eHlpq3hPlsnNft3wvsdQLNw3e0Jr05J6zwIjCf4Anw4adkRwOHhpaMtZrYFuBH4RopjcATwvaS2Q4B2qT5rMndfBlxDkMTXm9ljZlbqOJQjeduJX+qJX+jJ/gY0Do/vEUBv4FkAM2tiZn82s1Xh/98coGV4RlaRjkBVk0HyZykzRoLj/l9Jx70jcHgajqlkCSUoKc96gstf3wT+N2H++cBZBGcoLYDO4XxLaFNytmRmzQguj32RYh+PADOBju7eArg7aTtlcvdmCa9/JS1+muDS2PLwzCPRamCFu7dMeOW5+7cobTXwUFLbpu5+W2IoFcT5iLsPIfjSdeBX4aKvgCYJTQ9LtXrS9JPA0PCe0XcpI0G5+37gCYJkdj7wQni2BPBfBJdMB7h7c+CkcH6U474aOCrF/K/Cf8v7PAd8lgpiXA1MTTruTdz90XDdso6p1CJKUFIud/+C4FLSKDP7fTg7j+ByyyaCL6T/SbHqt8xsiJk1JLgX9a67pzrTyAO+dPciM+tP8EWVjri/CuNO1V38PWBbeKO9cXhGeJyZ9UvR9mHg22Z2WtguN+zMELVTwTFmdoqZNQKKCO737AsXLyQ4ToeY2WEEZwUVfa4NwJsE95BWuPuScpo/ApxLcD8wMZHlhXFsMbNDgJuifJbQvQQdHI62wPFmdmgY1xrgB+FxmkDqRBY1xnuAy8KzKzOzphZ0qMmr4JhKLaIEJRUKE8spwBgz+yXB5bNVBF9IiwnuUyV7hOCL70vgBIIvoFR+DEwxs0JgMsFf1OmKe36qexPuvg/4NsElpRXARoIv3hYp2q4mOFu8EdhA8Jf99UT/3WkE3Bbu498EHQRuDJc9RHBPbSXwKvB4xG0+QnD2WtblveLY3yU4szmcoFNLsT8AjcOY5gEvR9wvwO8I/o9eBbYB94XbguC+0fUEf7j0AN5OtYEoMYadci4h6LyxGVhGcMkWyj+mUouYChZKupnZdKDA3SdlOhYRyV46gxIRkVhSghIRkVjSJT4REYklnUGJiEgsZWzAxtatW3vnzp0ztXsREakhCxYs2OjubSpueaCMJajOnTszf37y8G4iIlLbmFnyw/KR6BKfiIjEkhKUiIjEkhKUiIjEUoX3oMzsfuAMYL27H5diuQF3EJQf2AGMd/cPktuJSN2wZ88eCgoKKCoqynQoUsNyc3Pp0KEDOTk5adlelE4S0wnGw0quqVNsNEHhuqMJCpL9KfxXROqggoIC8vLy6Ny5MwlVWqSWc3c2bdpEQUEBXbp0Scs2K0xQ7j7HzDqX0+Qs4EEPnvidZ2Ytzaydu68tb7s7137Gov8ZEjnQuY2H8UaTVNUQEgLp3Z7zBySXBRKRmlRUVKTkVAeZGYceeigbNmxI2zbTcQ+qPQcWIisI55ViZpea2Xwzm1+ZESw671nO4J2zy22zeO02/rpwTeRtikj1UXKqm9L9/56O56BSRZQy+7j7NGAaQH5+vve48f9F28MDp9MDePyHJ5bZ5Nw/vxNtWyIikhXScQZVQEL1VKADqSunVrstO/cwc+GaCl+zl67PRHgiUkMmTJhA27ZtOe64Uv26yvTmm29iZtx3330l8z788EPMjN/85jflrnv33Xfz4INl3aYPLFy4kFmzZlUYx/z587nqqquiBV2B6dOnc8UVV6RlW5mQjgQ1E7ggrHo5ENha0f2n6rJv337a5OVW+CrcuScT4YlIDRk/fjwvv1yZOoyBnj178vjjX9eNfOyxx+jVq1eF61122WVccMEF5baJmqDy8/O58847Kw62DqgwQZnZo8A7wDFmVmBmF5nZZWZ2WdhkFrCcoOLlPQQVUkVEMuakk07ikEMOqfR6nTp1oqioiHXr1uHuvPzyy4wePbpk+T333EO/fv3o1asX55xzDjt27ADg5ptvLjnLGjp0KD/72c/o378/Xbt25a233mL37t1MnjyZxx9/nN69e/P444/z3nvvMWjQIPr06cOgQYNYunQpEJzJnXHGGSXbnTBhAkOHDuXII488IHE9/PDD9O/fn969e/OjH/2IffuCqvcPPPAAXbt25eSTT2bu3LlVO4AxEaUX39gKljvwf9IWkYjUGrc8v4jFX2xL6za7H96cm77do9Lr3X777cyYMaPU/JNOOumAL/4xY8bw5JNP0qdPH/r27UujRo1Klp199tlccsklAEyaNIn77ruPK6+8stQ29+7dy3vvvcesWbO45ZZbeP3115kyZQrz58/nj3/8IwDbtm1jzpw5NGjQgNdff50bb7yRp59+utS2PvvsM2bPnk1hYSHHHHMMl19+OcuWLePxxx9n7ty55OTk8OMf/5gZM2YwYsQIbrrpJhYsWECLFi0YNmwYffr0qfSxiouMDRYrIlKTrr/+eq6//voK233/+9/n3HPP5bPPPmPs2LG8/fbbJcs+/fRTJk2axJYtW9i+fTunnXZaym2cffbZAJxwwgmsXLkyZZutW7dy4YUX8s9//hMzY8+e1LceTj/9dBo1akSjRo1o27Yt69at44033mDBggX069cPgJ07d9K2bVveffddhg4dSps2wcDh5557Lv/4xz8q/MxxlT0JaueX8MlTqZfltgCa1Wg4IlKxqpzpVJeoZ1CHHXYYOTk5vPbaa9xxxx0HJKjx48fz3HPP0atXL6ZPn86bb76Zcl/FZ13169dn7969Kdv893//N8OGDePZZ59l5cqVDB06tNxtJW7P3bnwwgv55S9/eUDb5557rlZ18c+eBLV/HzT7Rupl29dRmQRVtHcfMyvxzFRe4xyGHdM2cnsRiZ+oZ1AAU6ZMYf369dSvX/+A+YWFhbRr1449e/YwY8YM2rdP+chnSnl5eRQWFpZMb926tWT96dOnR94OwPDhwznrrLO49tpradu2LV9++SWFhYUMGDCAq6++mk2bNtG8eXOefPLJSJ084qpODhbbsVXTSL391OtPJDuNHTuWE088kaVLl9KhQ4cDuo5HMWjQIL7zne+Umn/rrbcyYMAARowYwbHHHlupbQ4bNozFixeXdJL46U9/ys9//nMGDx5c0sEhqu7du/OLX/yCkSNHcvzxxzNixAjWrl1Lu3btuPnmmznxxBM59dRT6du3b6W2GzdWmREd0ik/P98jFyx84HT4agOc/tvUy7ev49y327Np+y5u/U7P9AUZ2lBYxJm9o/+lJFKXLVmyhG7dumU6DMmQVP//ZrbA3fMru606eQYlIiLxpwQlIiKxpAQlIiKxlD29+CJYs2UnU15YVGG7wUe1Zni3MnoEiohILNSaM6izerenfcvGFbZbtWkHcz/fWAMRiYjIwag1Z1DnD+hEs0b1aZOXW267KGdYySrz3JSemRIRSY9acwZVnSrz3JSemRLJrNWrVzNs2DC6detGjx49uOOOOyKvW1xy4/nnny+Zd8YZZ5Q5YkSxP/zhDyUDx6bboEGDKmwTdf8XX3wxixcvTkdYdO7cmY0bq/dqlBKUiNQqDRo04Le//S1Llixh3rx53HXXXZX6Uu7QoQNTp06t1D6rM0ElDrV0sPu/99576d69ezrCqhFKUCJSq7Rr165kBIW8vDy6devGmjXRhzbr1asXLVq04LXXXiu17I033qBPnz707NmTCRMmsGvXLu68806++OILhg0bxrBhw0qt07lzZ2688UZOPPFE8vPz+eCDDzjttNM46qijuPvuuwHYvn07w4cPp2/fvvTs2ZO//vWvJes3axYM4/bmm28ydOhQxowZw7HHHsu4ceNw95T7v/zyy8nPz6dHjx7cdNNNJdsaOnQoxQMkNGvWjIkTJ9KrVy8GDhzIunXrANiwYQPnnHMO/fr1o1+/fiUlOzZt2sTIkSPp06cPP/rRj6iJQR5qzT0oEYmhl26Af3+S3m0e1hNG3xap6cqVK/nwww8ZMGAAEH3A2EmTJjFp0iRGjBhRMq+oqIjx48fzxhtv0LVrVy644AL+9Kc/cc011/C73/2O2bNn07p165RxdOzYkXfeeYdrr72W8ePHM3fuXIqKiujRoweXXXYZubm5PPvsszRv3pyNGzcycOBAzjzzzFIDv3744YcsWrSIww8/nMGDBzN37lyuuuqqUvufOnUqhxxyCPv27WP48OF8/PHHHH/88Qds66uvvmLgwIFMnTqVn/70p9xzzz1MmjSJq6++mmuvvZYhQ4bwr3/9i9NOO40lS5Zwyy23MGTIECZPnsyLL77ItGnTIv0fHAwlKBGplbZv384555zDH/7wB5o3bw5EHzD2m9/8JgBvvfVWybylS5fSpUsXunbtCsCFF17IXXfdxTXXXFPh9s4880wgqNi7fft28vLyyMvLIzc3ly1bttC0aVNuvPFG5syZQ7169VizZg3r1q3jsMMOO2A7/fv3p0OHDgD07t2blStXMmTIkFL7e+KJJ5g2bRp79+5l7dq1LF68uFSCatiwYUlhxBNOOKHkjPH1118/4JLotm3bKCwsZM6cOTzzzDNAUAKkVatWFX7ug1U7EtSeIvjkKQ5fvZnmjXNSNtnbsDlb259cw4GJ1HERz3TSbc+ePZxzzjmMGzeupDYTRD+DApg4cSJTp06lQYPga/JgLmkVl8yoV6/eAeUz6tWrx969e5kxYwYbNmxgwYIF5OTk0LlzZ4qKisrcDpRdymPFihX85je/4f3336dVq1aMHz8+5bZycnJKztASt7V//37eeecdGjcu/dhOTZfyqB33oFodAc2+wa7c1uxp3Cblq8Hur6t6bijclcFgRaQ6uTsXXXQR3bp147rrrjtg2fXXX8/ChQtLvZKTE8DIkSPZvHkzH330EQDHHnssK1euZNmyZQA89NBDnHxy8EdvcimNytq6dStt27YlJyeH2bNns2rVqkqtn7j/bdu20bRpU1q0aMG6det46aWXKrWtkSNHllT9BVi4cCEQJPHi5P7SSy+xefPmSm23KmpHgqqkjdt3ZzoEEakmc+fO5aGHHuJvf/sbvXv3pnfv3syaNatK25o4cSIFBQUA5Obm8sADD/C9732Pnj17Uq9ePS677DIALr30UkaPHp2yk0QU48aNY/78+eTn5zNjxoxKl/JI3H+vXr3o06cPPXr0YMKECQwePLhS27rzzjuZP38+xx9/PN27dy/pyHHTTTcxZ84c+vbty6uvvkqnTp0qtd2qyJ5yG18sgEOOSr38yKHQdRRzl22kZZOGKZvk7NzApi7fZsoLi1iytpBHLxlYpbgrotIcUtep3Ebdls5yG9lxD6rnmKAeVCpfrgj+7Tqq5uIpR2Wr9ZZHo1KISF2WHQkq/4fQKC91yfeXb6j5eMrRsVXTtG1rQ2HpG5siInVFpARlZqOAO4D6wL3uflvS8k7AX4CWYZsb3L1qF31rSGXG5NPo5yKV4+413uNLMi/dt4wqTFBmVh+4CxgBFADvm9lMd08cO2QS8IS7/8nMugOzgM5pjbQKZiyFR/5RPNUGmFeybMna0j1uWjdrSJu8RgfMW7VpB7BRCUokotzcXDZt2sShhx6qJFWHuDubNm0iN7f8AbsrI8oZVH9gmbsvBzCzx4CzgMQE5UDz8H0L4Iu0RXgQxh0TvODrThIAY++ZF7mTRFVGPxepyzp06EBBQQEbNpRx31hqrdzc3JIHidMhSoJqD6xOmC4ABiS1uRl41cyuBJoCp6YlOhHJOjk5OXTp0iXTYUgtECVBpTpHT77QOBaY7u6/NbMTgYfM7Dh333/AhswuBS4FaqQPfaJ6+4o4dEXxEPptEt6XL2dnC4DI7TVihYhIekRJUAVAx4TpDpS+hHcRMArA3d8xs1ygNbA+sZG7TwOmQfAcVBVjrpJdzToeML2ncZtI63n9yrXP2anLGiIi6RBlJIn3gaPNrIuZNQTOA2YmtfkXMBzAzLoBuYC+qUVEpMoqTFDuvhe4AngFWELQW2+RmU0xszPDZv8FXGJmHwGPAuM9U0NURHB+10xHICIiFYn0HFT4TNOspHmTE94vBio34FMGFffsExGR+KqTg8WKiEj8KUGJiEgsZcdYfHVUeQPPaiBZEantlKBirLyBZzWQrIjUdrrEJyIisaQEJSIisaRLfGl24JBKZdOQSCIi5VOCimD5Vrjh7WhtT27fkdFHVNxOQyKJiJRPCaoCJ7eP3nb51uDfKAlKRETKpwRVgdFHRE84Uc+yRESkYuokISIisaQEJSIisaQEJSIisaQEJSIisaQEJSIisVQ7evF9uQJevoHjdu6hQf3UOXfrYYPY3GF4DQcmIiJVlf0J6sihFTbJLVwFoAQlIpJFsj9BdR0VvIBPl22kZZOGpZp0nn9rTUdV7corxVFTVPJDRKpT9ieoLFU8Zt+DnzfhgqN2HLAsyjh95ZXiqCkq+SEi1UkJKkN2NesIwMMrYOxxByYbjdMnIqJefCIiElNKUCIiEktKUCIiEkuR7kGZ2SjgDqA+cK+735aizfeBmwEHPnL389MYZ41p8/lTbDhqTJXXr0ztqGLJ7W1fC/YsWsTgo1ozvNs3qhyLiEg2qzBBmVl94C5gBFAAvG9mM919cUKbo4GfA4PdfbOZZW3f47bLn6lygiqvdtS6HbB+Z+pln2xKntOQnMLtAEpQIlJnRTmD6g8sc/flAGb2GHAWsDihzSXAXe6+GcDd16c70GxQmdpRxU5/Hl789oHzcnZu4OpFR6YvMBGRLBTlHlR7YHXCdEE4L1FXoKuZzTWzeeElwVLM7FIzm29m8zdsUFdqEREpW5QEZSnmedJ0A+BoYCgwFrjXzFqWWsl9mrvnu3t+mzZtKhuriIjUIVEu8RUAHROmOwBfpGgzz933ACvMbClBwno/LVGmQW7hqshDHh3M0EgalFZEJD2iJKj3gaPNrAuwBjgPSO6h9xzBmdN0M2tNcMlveToDPRhbDxtUal7Ozg00LNqYsn3TzUtKzdud25o9jcs/69OgtCIi6VNhgnL3vWZ2BfAKQTfz+919kZlNAea7+8xw2UgzWwzsA65391J90zJlc4fhkZNGj9fOZ9GIR6q0n6qceZ3ftUq7EhGp9SI9B+Xus4BZSfMmJ7x34LrwJZUw7phMRyAiEk8aLDaG6u0rImdncAJ66IrnD1gWZaRzEZHaQAkqhnY164jXD94n3/fSSOciUldoLD4REYklJagk6488O9MhiIgISlClHMxAsSIikj5KUCIiEkvqJCFVVrR3HzMXrsl0GCISc/UaNWtelfWUoKTKOrZqmukQRCQb1KtXv0qrpTsOERGRdKhVZ1BNGjVgy47dFbbbvW8/bfNyayAiERGpqlqVoPp0LFXhI6W5y1IPEhs3qcrHqxy8iNQVtSpBxUHUsh4VleUor3z8qk07gI1KUCJSqylBpVGqsh6pRCnLUVb5+JydW7l60aFVik9EJJtkT4LKbQHb16VetqcIWqX4Nq9hUct6HExBRBGRuiJ7EtTRI8pe9slTNReHiIjUCHUzFxGRWMqeMygByq8VlQ6qNyUicaEElWXKqxWVDqo3JSJxUScTVNQHeovpwV4RkZpXJxNU1Ad6i2XLg70iIrWJOkmIiEgsKUGJiEgsRUpQZjbKzJaa2TIzu6GcdmPMzM0sP30hiohIXVRhgjKz+sBdwGigOzDWzLqnaJcHXAW8m+4gRUSk7olyBtUfWObuy919N/AYcFaKdrcCvwaK0hifiIjUUVESVHtgdcJ0QTivhJn1ATq6+wtpjE1i4KkFqytuJCJSDaJ0M7cU87xkoVk94PfA+Ao3ZHYpcClAp06dokVYS5VXlqOiUhw16ekP1jDmhI6ZDkNE6qAoZ1AFQOI3VAfgi4TpPOA44E0zWwkMBGam6ijh7tPcPd/d89u0Sf8oCNli62GDKMpLPfp6buEqWvz77ZTLRETqkihnUO8DR5tZF2ANcB5wfvFCd98KtC6eNrM3gZ+4+/z0hlp7lFeWQ6U4REQCFZ5Bufte4ArgFWAJ8IS7LzKzKWZ2ZnUHKCIidVOkoY7cfRYwK2ne5DLaDj34sKQiy7fCDdVwJdD2tWDPokUHzJvywqKUbQcf1Vpl50Wk2tTJsfiy3cntK25TWet2wPqdAA1hS+EBy5asLSzVPqd+0HdGCUpEqosSVBYafUTwqg45Ozewqcu3S6bH3jOPRy8ZWKpdWWdVIiLpogQVQWXKc6g0h4hIeihBRVCZ8hzZXpqj3r6ipEq9dfdxABHJrNqRoHJbwPZ15bfZUwStqum6WC2yq5keyhWReKgdCeroERW3+eSp6o9DRETSRvWgpFw/6PJVpkMQkTpKCUrKdcFROzIdgojUUUpQIiISS0pQIiISS0pQIiISS7WjF18tky21okREqpMSVMxsPWxQmctyC1cBKEGJSJ2gBBUzqhUlIhLQPSgREYmlunMGFWU4pDRo8tVaaHJkte9HRKS2qzsJKspwSGnQ8Iv7I498XhGNjC4idVndSVA1pMfhzaFZ67RsK9tHRhcRORi6ByUiIrGkBCUiIrGkBCUiIrGkBCUiIrGkBCUiIrEUKUGZ2SgzW2pmy8zshhTLrzOzxWb2sZm9YWaqrS4iIgDUy23WqkrrVdTAzOoDdwGjge7AWDPrntTsQyDf3Y8HngJ+XZVgRESk9qnXqEn1JCigP7DM3Ze7+27gMeCsxAbuPtvdi0uvzgM6VCUYERGRYlEe1G0PrE6YLgAGlNP+IuClVAvM7FLgUoBOnTpFDFESlVeKIx1UzkNE4iJKgrIU8zxlQ7MfAPnAyamWu/s0YBpAfn5+ym1kvahj/u0pglaVu1VXXimOdFA5DxGJkygJqgDomDDdAfgiuZGZnQpMBE52913pCS8LRR3z75OnKr3p8kpxpENlz8xWbdrBlBcWVVM0IlLXRUlQ7wNHm1kXYA1wHnB+YgMz6wP8GRjl7uvTHqXEzuCjWgMaK1BEvrahcBcbt6dnsGyIkKDcfa+ZXQG8AtQH7nf3RWY2BZjv7jOB24FmwJNmBvAvdz8zbVFK7Azv9g2Gd/tGpsMQkSww9C9VWy/SaObuPguYlTRvcsL7U6u2eylPk0YNyizdoVIcIlLbqdxGjPXp2LLMZSrFISK1nYY6EhGRWFKCEhGRWFKCEhGRarV/147NVVlPCUpERKrV/qLtSlAiIlJ7qBdfppQ3JFIVhkESEaltlKAypbwhkaowDJKISG2jS3wiIhJLOoOSAySX89h+6PFs6vLtDEYkInWVEpSUSC7nkVu4Ctu3J0PRiEhdpwSVpcobp6+qthzyTTjkmyXTx33yS+rt38OGwqKU7Yv27qNjq6ZpjUFEpJgSVJYqb5y+tFmWA/vgzN7tUy6euXBN9ccgInWWOkmIiEgs6QwqjvSMlIiIElQs6RkpERFd4hMRkXhSghIRkVhSghIRkVhSghIRkVhSghIRkVhSL75sU14X9HRQN3YRiYlICcrMRgF3APWBe939tqTljYAHgROATcC57r4yvaEKUH4X9HRQN3YRiYkKL/GZWX3gLmA00B0Ya2bdk5pdBGx29/8Afg/8Kt2BiohI3RLlDKo/sMzdlwOY2WPAWcDihDZnATeH758C/g68Ja4AAA93SURBVGhm5u6exlglZvIa55Q5kKyISIn9+/dVZbUoCao9sDphugAYUFYbd99rZluBQ4GNiY3M7FLgUoBOnTpVJV6pbon3uPIOg/qNymw67Ji2NRSUiGSz/bu2b6vKelESlKWYl3xmFKUN7j4NmAaQn5+vs6s4SrzH1XNM5uIQkTovSjfzAqBjwnQH4Iuy2phZA6AF8GU6AhQRkbopSoJ6HzjazLqYWUPgPGBmUpuZwIXh+zHA33T/SUREDkaFl/jCe0pXAK8QdDO/390XmdkUYL67zwTuAx4ys2UEZ07nVWfQIiJS+0V6DsrdZwGzkuZNTnhfBHwvvaGJiEhdpqGOREQklpSgREQklixTfRnMrBBYmpGdH7zWJD3jlWWyOX7FnhmKPXOyOf7i2I9w9zaVXTmTg8Uudff8DO6/ysxsfrbGDtkdv2LPDMWeOdkc/8HGrkt8IiISS0pQIiISS5lMUNMyuO+Dlc2xQ3bHr9gzQ7FnTjbHf1CxZ6yThIiISHl0iU9ERGJJCUpERGIpIwnKzEaZ2VIzW2ZmN2QihqjMrKOZzTazJWa2yMyuDucfYmavmdk/w39bZTrWsphZfTP70MxeCKe7mNm7YeyPh4MAx46ZtTSzp8zss/D4n5gtx93Mrg1/Xj41s0fNLDfOx93M7jez9Wb2acK8lMfaAneGv78fm1nfzEVeZuy3hz83H5vZs2bWMmHZz8PYl5rZaZmJuiSWUrEnLPuJmbmZtQ6nY3/cw/lXhsd2kZn9OmF+5Y+7u9foi2DA2c+BI4GGwEdA95qOoxLxtgP6hu/zgH8A3YFfAzeE828AfpXpWMv5DNcBjwAvhNNPAOeF7+8GLs90jGXE/Rfg4vB9Q6BlNhx3ggKeK4DGCcd7fJyPO3AS0Bf4NGFeymMNfAt4iaAO3EDg3RjGPhJoEL7/VULs3cPvnEZAl/C7qH6cYg/ndyQYoHsV0DqLjvsw4HWgUTjd9mCOeybOoEpKyLv7bqC4hHwsuftad/8gfF8ILCH4AjqL4AuU8N/vZCbC8plZB+B04N5w2oBTgKfCJrGM3cyaE/wC3Afg7rvdfQtZctwJHoJvHNZHawKsJcbH3d3nULqGW1nH+izgQQ/MA1qaWbuaibS0VLG7+6vuvjecnEdQxw6C2B9z913uvgJYRvCdlBFlHHeA3wM/5cDCr7E/7sDlwG3uvitssz6cX6XjnokElaqEfPsMxFFpZtYZ6AO8C3zD3ddCkMSAuNY//wPBD/r+cPpQYEvCL29cj/+RwAbggfDy5L1m1pQsOO7uvgb4DfAvgsS0FVhAdhz3RGUd62z7HZ5AcOYBWRC7mZ0JrHH3j5IWxT52oCvwzfBS9t/NrF84v0qxZyJBRSoPHzdm1gx4GrjG3bdlOp4ozOwMYL27L0icnaJpHI9/A4LLB39y9z7AVwSXmWIvvFdzFsGljMOBpsDoFE3jeNyjyJafIcxsIrAXmFE8K0Wz2MRuZk2AicDkVItTzItN7KEGQCuCS5DXA0+EV22qFHsmElSUEvKxYmY5BMlphrs/E85eV3x6Hf67vqz1M2gwcKaZrSS4lHoKwRlVy/DSE8T3+BcABe7+bjj9FEHCyobjfiqwwt03uPse4BlgENlx3BOVdayz4nfYzC4EzgDGeXgjhPjHfhTBHzYfhb+3HYAPzOww4h87BDE+E16GfI/gyk1rqhh7JhJUlBLysRFm//uAJe7+u4RFiWXuLwT+WtOxVcTdf+7uHdy9M8Fx/pu7jwNmA2PCZnGN/d/AajM7Jpw1HFhMFhx3gkt7A82sSfjzUxx77I97krKO9UzggrBX2UBga/GlwLgws1HAz4Az3X1HwqKZwHlm1sjMugBHA+9lIsZU3P0Td2/r7p3D39sCgk5a/yYLjjvwHMEfwphZV4LOTRup6nHPUO+PbxH0hvscmJiJGCoR6xCCU9GPgYXh61sE93LeAP4Z/ntIpmOt4HMM5etefEeGPxzLgCcJe9zE7QX0BuaHx/45gksHWXHcgVuAz4BPgYcIei/F9rgDjxLcL9tD8KV4UVnHmuByzV3h7+8nQH4MY19GcM+j+Hf27oT2E8PYlwKj4xZ70vKVfN2LLxuOe0Pg4fDn/gPglIM57hrqSEREYkkjSYiISCwpQYmISCwpQYmISCwpQYmISCwpQYmISCwpQYmISCwpQUmtYUF5jh8nTB9uZk+Vt04ltn2zmf0kfD/FzE5N03bbWVgGpTqY2fZKtH3dYlq+ROomJSipTVoCJQnK3b9w9zHltK8Sd5/s7q+naXPXAfekaVsH6yESjp9IpilBSW1yG3CUmS0MC9Z1Li6mZmbjzew5M3vezFaY2RVmdl04Uvo8MzskbHeUmb1sZgvM7C0zOzZ5J2Y23czGhO9XmtktZvaBmX1S3N7MmoYF3d4P91FWSZlzgJfDdWaZ2fHh+w/NbHL4/lYzuzh8f324zY/N7JaEmH5gZu+Fn/3PZlY/KebWZvaOmZ0enrXNCdt+ambfDJvNBMZW8diLpJ0SlNQmNwCfu3tvd78+xfLjgPMJ6tBMBXZ4MFL6O8AFYZtpwJXufgLwE+B/I+x3o7v3Bf4UrgPBsC5/c/d+BEXcbg/LhZQIxyTb7GHtHGAOQamC5gQjcA8O5w8B3jKzkQRjmPUnGAbqBDM7ycy6AecCg929N7APGJewn28ALwKT3f3F8Bi8ErbtRTAUEO6+GWhkZodG+Mwi1a5BxU1Eao3ZHhSdLDSzrcDz4fxPgOPDkiqDgCeDMV6BYAy9ihSPcL8AODt8P5JgJPnihJULdCIoeFmsHUHNq2JvAVcRVON9ERgRll/o7O5LzeyScLsfhu2bESSs44ETgPfDuBvz9cjjOQTj6P0fd/97OO994P5wlP7n3H1hQgzrCUqEbIrwuUWqlRKU1CW7Et7vT5jeT/C7UI+gqGDvKm53H1//ThlwjrsvLWe9nQSJq9j7QD6wHHiNoEzBJQSJr3ibv3T3PyduxMyuBP7i7j9PsY+94fqnAX+HoBKqmZ1EUGn5ITO73d0fDNvnhnGJZJwu8UltUgjkVXVlDwpRrjCz70FQasXMelVxc68AV4blNjCzPina/APonLD/3QQjcH+foEz5WwSXDN9K2OaE8EwPM2tvZm0JzpDGhO8xs0PM7IjizRJUlD3WzG4Ilx9BUMjyHoJSMn2LPy9wGMEI2iIZpwQltYa7bwLmhjf+b6/iZsYBF5nZR8Aigsq4VXErweW1j8OOGremiPcr4HMz+4+E2W8B6zyoYfQWQWG3t8L2rwKPAO+Y2ScERRzz3H0xMAl41cw+Jjj7apewn30E9cCGhd3whwILzexDgk4ad4RNTwDm+ddl6UUySuU2RDLIzL4LnODuk2IQyx3ATHd/I9OxiEAazqDM7Ltm5ondccPuvedXYhtvV7A88sOGItnE3Z8lPpfUPlVykjg56DMoM3uC4HLCG+5+czhvKPATdz+jgnXrh5cfKtrHdndvdlCBiohIVjmoM6jwZu1gglK/5yUsuo3geY6FZnZt0jpDzWy2mT1C0L235AypnAcIi9ctedjwYOIWEZH4O9hu5t8BXnb3f5jZl2bW190/IHhgsrwzqP7Ace6+Iml+8QOEU8Mn4ZsULwgfNpwJTHL31w4ybhERibmDvQc1FngsfP8Y0YdJeS9FcoLgOZAfmtnNQM/woUr4+mHDnyo5iYjUDVVOUOFwKKcA95rZSuB64Nzi5z4q8FWqme4+BzgJWEPwAGHx8DOJDxuKiEgdcDBnUGOAB939CHfv7O4dCYZoGUIVH5gs6wFCUjxsKCIitdvBJKixwLNJ854muI/0MbDXzD5K7iRRgaGkfoAw1cOGIiJSi+lBXRERiSUNdSQiIrGkBCUiIrGkBCUiIrGkBCUiIrGkBCUiIrGkBCUiIrGUtgRlZveb2fqwOFvxvO+Z2SIz229m+enaV00ys2PCwWuLX9vM7JpMxxWFmeWa2Xvh82iLzOyWTMcUVaqfp2ySzfEr9szI5tiheuJP5xnUdGBU0rxPgbOBOWncT41y96Xu3tvdexNUHN1B6QeU42oXcIq79wJ6A6PMbGCGY4pqOqV/nrLJdLI3/uko9kyYTvbGDtUQf9oSVDiO3pdJ85a4+9J07SMGhgOfu/uqTAcShQeKiz3mhK+seDI71c9TNsnm+BV7ZmRz7FA98eseVOWcBzya6SAqw8zqm9lCYD3wmru/m+mYRESiUIKKyMwaAmcCT2Y6lspw933h5ckOQH8zOy7TMYmIRKEEFd1o4AN3X5fpQKrC3bcAb5Ld17hFpA5RgopuLNl3ea+NmbUM3zcGTgU+y2xUIiLRpLOb+aPAO8AxZlZgZheZ2XfNrAA4EXjRzF5J1/5qkpk1AUYAz2Q6lkpqB8w2s48JqhW/5u4vZDimSFL9PGU6psrI5vgVe2Zkc+xQPfGr3IaIiMSSLvGJiEgsKUGJiEgsKUGJiEgsVfdYfLeb2Wdm9rGZPVvcoyzbmFlLM3sq/CxLzOzETMcUhZl1NLPZYcyLzOzqTMdUGWY2ysyWmtkyM7sh0/FURjbHDtkdv2LPjGqJ3d3T8gJOAvoCnybMGwk0CN//CvhVuvZXky/gL8DF4fuGQMtMxxQx7nZA3/B9HvAPoHum44oYe33gc+DI8Jh/pNgVv2KP56u6Yq/usfhedfe94eQ8gtEMsoqZNSdIvvcBuPtuDx56jT13X+vuH4TvC4ElQPvMRhVZf2CZuy93993AY8BZGY4pqmyOHbI7fsWeGdUSe03eg5oAvFSD+0uXI4ENwANm9qGZ3WtmTTMdVGWZWWegD5AtY/G1B1YnTBeQPck1m2OH7I5fsWdGtcReIwnKzCYCe4EZNbG/NGtAcOnyT+7eB/gKyLZrw82Ap4Fr3H1bpuOJyFLMy5aH9rI5dsju+BV7ZlRL7NWeoMzsQuAMYJyHFyuzTAFQ4F+PAv4UQcLKCmaWQ5CcZrh7No2EUQB0TJjuAHyRoVgqK5tjh+yOX7FnRrXEXq0JysxGAT8DznT3HdW5r+ri7v8GVpvZMeGs4cDiDIYUmZkZwb2zJe7+u0zHU0nvA0ebWZdwJPnzgJkZjimqbI4dsjt+xZ4Z1RJ7g4MOKxSOwzQUaB2Ov3cT8HOgEfBa8F3JPHe/LF37rEFXAjPCA78c+GGG44lqMPCfwCdhTSiAG919VgZjisTd95rZFcArBD2E7nf3RRkOK5Jsjh2yO37FnhnVFbvG4hMRkVjSSBIiIhJLSlAiIhJLSlAiIhJLSlAiIhJLSlAiIhJLSlAiIhJLSlAiIhJL/x/Qg6xPjjeMFgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T1, T2 = df['time']  [df['group'] == 1], df['time']  [df['group'] == 2]\n",
    "E1, E2 = df['status'][df['group'] == 1], df['status'][df['group'] == 2]\n",
    "kmf1, kmf2 = fit_plot(T1, T2, E1, E2, 'Kaplan-Meier survival curves', 'weeks', '1=Maintained', '2=Not maintained')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vertical ticks on the graphs indicate censored events. At the bottom, the number of patients \"at risk\" is shown. The study concerned 11 patients in group 1 with one patient censored after 161 weeks and 12 patients in group 2. A Kaplan-Meier curve estimates the survival probability as a function of time (duration). \n",
    "\n",
    "In this example, the two curves do not appear to differ very much. To analyze this more precisely one performs a [logrank test](https://en.wikipedia.org/wiki/Logrank_test):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06533932204050484"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lifelines.statistics.logrank_test(T1, T2, E1, E2).p_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null hypothesis is that survival is the same for both groups. Since $p=0.065$ is not particularly small (e.g., not below $\\alpha=0.05$), the null hypothesis is not strongly rejected. In other words, the logrank test also says that the curves do not differ significantly---with the usual caveat of small sample sizes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Privacy-Preserving Survival Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a multiparty setting, each party holds a private dataset and the goal is to perform a survival analysis on the *union* of all the datasets. To obtain the secure union of the datasets in an efficient way, each private dataset will be represented using two survival tables, one survival table per group of patients. \n",
    "\n",
    "The survival tables for the aml dataset are available from the Kaplan-Meier fitters `kmf1` and `kmf2` output by the call to `fit_plot()` above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>removed</th>\n",
       "      <th>observed</th>\n",
       "      <th>censored</th>\n",
       "      <th>entrance</th>\n",
       "      <th>at_risk</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>event_at</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>9.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>13.0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>18.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>23.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>28.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>31.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>34.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>45.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>48.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>161.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          removed  observed  censored  entrance  at_risk\n",
       "event_at                                                \n",
       "0.0             0         0         0        11       11\n",
       "9.0             1         1         0         0       11\n",
       "13.0            2         1         1         0       10\n",
       "18.0            1         1         0         0        8\n",
       "23.0            1         1         0         0        7\n",
       "28.0            1         0         1         0        6\n",
       "31.0            1         1         0         0        5\n",
       "34.0            1         1         0         0        4\n",
       "45.0            1         0         1         0        3\n",
       "48.0            1         1         0         0        2\n",
       "161.0           1         0         1         0        1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>removed</th>\n",
       "      <th>observed</th>\n",
       "      <th>censored</th>\n",
       "      <th>entrance</th>\n",
       "      <th>at_risk</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>event_at</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>5.0</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>8.0</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>12.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>16.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>23.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>27.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>30.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>33.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>43.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>45.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          removed  observed  censored  entrance  at_risk\n",
       "event_at                                                \n",
       "0.0             0         0         0        12       12\n",
       "5.0             2         2         0         0       12\n",
       "8.0             2         2         0         0       10\n",
       "12.0            1         1         0         0        8\n",
       "16.0            1         0         1         0        7\n",
       "23.0            1         1         0         0        6\n",
       "27.0            1         1         0         0        5\n",
       "30.0            1         1         0         0        4\n",
       "33.0            1         1         0         0        3\n",
       "43.0            1         1         0         0        2\n",
       "45.0            1         1         0         0        1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(kmf1.event_table, kmf2.event_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `mpc.input()`, each party will secret-share its pair of survival tables with all the other parties. To allow for a simple union of the survival tables, however, we modify the representation of the survival tables as follows. \n",
    "\n",
    "First, we determine when the last event happens, which is at $t=161$ in the example. We compute $maxT$ securely as the maximum over all parties, using the following MPyC one-liner:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "161"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "secfxp = mpc.SecFxp(64)  # logrank tests will require fixed-point arithmetic\n",
    "maxT = int(await mpc.output(mpc.max(mpc.input(secfxp(int(df['time'].max()))))))\n",
    "timeline = range(1, maxT+1)\n",
    "max(timeline)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All parties use $1..maxT$ as the global timeline. Subsequently, each party pads its own survival tables to cover the entire timeline. This is accomplished by two calls to the function `events_to_table()`, which only keeps the essential information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>d1</th>\n",
       "      <th>n1</th>\n",
       "      <th>d2</th>\n",
       "      <th>n2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>2</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>2</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    d1  n1  d2  n2\n",
       "1    0  11   0  12\n",
       "2    0  11   0  12\n",
       "3    0  11   0  12\n",
       "4    0  11   0  12\n",
       "5    0  11   2  12\n",
       "6    0  11   0  10\n",
       "7    0  11   0  10\n",
       "8    0  11   2  10\n",
       "9    1  11   0   8\n",
       "10   0  10   0   8"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d1, n1 = events_to_table(maxT, T1, E1)\n",
    "d2, n2 = events_to_table(maxT, T2, E2)\n",
    "pd.DataFrame({'d1': d1, 'n1': n1, 'd2': d2, 'n2': n2}, index=timeline).head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Column `d1` records the number of observed events (\"deaths\") for group 1 on the entire timeline, and column `n1` records the number of patients \"at risk\". Similarly, for group 2. \n",
    "\n",
    "To obtain the secure union of the privately held datasets, we let all parties secret-share their survival tables and simply add all of these elementwise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "d1, n1, d2, n2 = (functools.reduce(mpc.vector_add, mpc.input(list(map(secfxp, _)))) for _ in (d1, n1, d2, n2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The joint dataset (union) is now represented by `d1, n1, d2, n2`. Note that these values are all secret-shared, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mpyc.sectypes.SecFxp64:32 at 0x1c355f048c8>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d1[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregate Kaplan-Meier Curves\n",
    "\n",
    "We now proceed to analyze the joint dataset in a privacy-preserving way by plotting *aggregated* versions of the Kaplan-Meier curves. The exact Kaplan-Meier curves would reveal too much information about the patients in the study. By aggregating the events over longer time intervals, the amount of information revealed by the curves is reduced. At the same time, however, the aggregated curves may still be helpful to see the overall results for the study---and in any case to serve as a sanity check.\n",
    "\n",
    "The function `aggregate()` securely adds all the observed (\"death\") events over intervals of a given length (stride). The aggregated values are all output publicly, and used to plot the curves via the function `events_from_table()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEOCAYAAADc94MzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwU1bn/8c+XdZDNBYgL6GAiKgRZHBGXKGhEvTGYRE003CtcjUZzjcYkGhcuosbExGz6i1cvasQFxSVqUHGPBi9xYVA0ImKIoowLAiqCimzP74+qGZuhZ6aBGapgvu/Xq1/TVXWq6pkz3f3MOVV9jiICMzOzvGmRdQBmZmbFOEGZmVkuOUGZmVkuOUGZmVkuOUGZmVkuOUGZmVkuOUFZbkkaK+nmDM47StL/bezzrgtJMyUNyTqOdSFpR0lLJbXcwOOUSwpJrRorNssnJ6ickvSEpA8ktc06lvUlaa6krzbRsYdIqipYbiPpLklTJXVqinM2hjTphqTTa63/Ubp+bCnHiYg+EfFEU8TYVCLizYjoEBGrso7FNg1OUDkkqRz4ChDA8CY8z2bxH2iaxO8CtgSGRcRHGYfUkFeBkbXWHZ+ub1JN+TffXF5PAEr48zFj/gPk0/HA08B4an2QSdpG0r2SPpI0TdLPC7ujJA2TNFvSYkn/I+lvkr6XbhuVtjB+L+l9YGy6/gRJs9IW20OSdirxeF+U9FdJiyQtlDRB0pbptpuAHYF7026ds9P1gyX9XdKHkl4o7KaS1DM9/hJJjwBdGqooSVsA9wKtga9FxMfp+kGSnkrP846kP0pqU7BfSDpd0mtp7JfV9YEk6XJJ89I6ny7pKwXbxkq6XdKNadwzJVU0EPY0YAtJfdJj9AHapesLz3uEpBnp7/B3SXsUbKtpnUpqIekcSf9K/xa3S9o63VbdHXaipDeBvxb5/bpIui89z/uSnqyui3TfLxWUHS/p5+nzIZKqJP1M0rvA9enr6IiC8q3S+h1YEEsrScdKqqwVx5mSJqXPvybp+bTO56nElmW6bw8lrekFaX38MV2/RpexanUVKum1uETSVOAT4LwGYmwr6TeS3pQ0X9LVkto1VKdWOldYPh0PTEgfh0r6QsG2K4GPgW1JkldNApPUBbgTOBfYBpgN7Fvr2HsDrwHdgEskfQM4D/gW0BV4Eri1xOMJ+CWwPbA70IM06UXEfwBvAl9Pu3V+LWkH4H7g58DWwE+BP0vqmh7vFmA6SWK6mLVbGbW1BR4AlgHDI+LTgm2rgDPTY+0DHAz8oNb+3wQqgIHAkcAJdZxnGtA/jfkW4A5JZQXbhwMTSVpwk4A/NhA3wE0kf2dIfs8bCzdKGgj8Cfg+Sd3/LzBJxbt8Twe+ARxI8rf4gOR1UuhAkr/RoUX2/wlQRfL3/wLJ66HUMdC2JamXnYCTSV47xxVsPxRYGBHP1dpvErCrpF0K1n2XpH4heY0fT1KnXwNOTV+r9VJyfes+4A2gHNiB5G9Tqv9If4+OwP9rIMZfAb1IXhtfSs81Jt22IXVq1SLCjxw9gP2BFUCXdPkV4Mz0ect0264F5X8O/F/6/HjgqYJtAuYB30uXRwFv1jrfA8CJBcstSP573Kmh4xWJ/RvA8wXLc4GvFiz/DLip1j4PkXxA7wisBNoXbLsFuLmOcw0hSUzLgaNKqNcfAXcXLAdwWMHyD4DHCurp/+o51gdAv/T5WODRgm29gU/r2XcscHP6+75J0vJ7kyS53wyMTctdBVxca9/ZwIG16xaYBRxcUG679HXSiuRDOoCd64npIuAvwJeKbIvC9SSt+p8X/A2WA2UF278ELAG2SJcnAGPS59WxtEqXby7YtkvhfkXi+APw+2LHqVVuH2BBHdvGFr6eisTzBHBRrX2KxkjyXvgY+GKtc7/eUJ36UfrDLaj8GQk8HBEL0+Vb+Lwl0ZXkQ2deQfnC59sXLkfyTqliTfNqLe8EXJ52RXwIvE/y5tuhoeNJ6iZpoqS3JH1E8maur1tuJ+CY6nOl59uf5AN1e+CDSLvoUm/UcyyAhcCxwA2S1mgZSOqVdrG8m8b2iyKxFdbFG2kMa5H0k7TranEac+dax3q34PknQFnajTVCSffmUkkPFB4zIt4E5qRx/TMiiv1dflKrrnrUEeNOwN0F5WaRtCALW961j1/osjSWh5V0eZ5TT9naFkTEsoLfa056/q8r6X4dzuctjtpu4fPW1neBeyLiEwBJe0t6PO2mWwycQgldviR19EZErFyH36FQ7XqqK8auJIlqekG9P5iuhw2rU0s5QeVI2n/9beDA9IP1XZJuqn6S+pH8Z7gS6F6wW4+C5+8UbpOkWmVh7W6GecD3I2LLgke7iPh7Ccf7ZXq8PSKiE/DvJMmtvnPdVOtc7SPi0vRcW0lqX1B+RxoQEXcBJwF3ShpasOkqktbnLmls59WKDdasux2Bt2sfX8n1pp+R/F22iogtgcVFjlUstgmRdG92iIjDixS5kaQr6MYi2+YBl9Sqqy0i4tY6yh5eq2xZRLxVGE49cS6JiJ9ExM7A14EfSzo43fwJyQdxtW1r717kkNXdfEcCL6dJq5iHgS6S+qflCxPZLSTdgD0iojNwNSXUOUld7KjiN2x8TP2/C6z9+9QV40LgU6BPQZ13jogO0GCdWomcoPLlGyT/+fYm6dfuT3Ld4Eng+Ehuz70LGCtpC0m78fl1DEiu7/SV9I30DfpfFH8TFroaOFefX7DvLOmYEo/XEVgKfJheXzqr1rHnAzsXLN9M8p/1oZJaSipTcqG9e0S8AVQCFyq5ZXx/kjd2g9IP7dOAv0jaryC2j4ClaT2dWmTXsyRtJakHcAZwW5EyHUn+KVgAtJI0Bmis29hvA4YBtxfZdg1wStqSkKT2Sm4c6Fik7NUk1xN3ApDUVdKRpQah5GaML6X/gHxE8hqsvhV8BvDd9O91GMm1rIZMTH+vU6m79UTayrmTpLWxNfBIweaOwPsRsUzSIJLWSymeJfln59K0zsoKXhMzgAOUfB+rM8m11XrVFWNErCb5G/1eUjcASTtUt+QbqFMrkRNUvowEro/k+yLvVj9ILrqPSJPEaSRdTO+SXGi/FfgMIO0WPAb4NbCIJNFVVm8vJiLuJrnYOzHtCnsJOLzE411IcoPBYpJkdletw/8SGJ12gfw07cY6kqQ1s4Dkv92z+Px1+F2SmzjeBy6geMuirt/jBpLWyP3pB9pP0+MtIfkgKZZ8/kJyU8aMNP7ripR5iOQ63ask3YDLqL+7rGQR8WlEPBpr3txRva2SpGX4R5JrXnNIro0VczlJa+NhSUtI7gDdex1C2QV4lOSfjaeA/4nPv2N1Bsk/Ch8CI4B7GjpYRLyTHmdfitd7oVuArwJ31OqW+wFwUfr7jKF4Ei927lVpvF8iubZXBXwn3fZIGs+LJH/3+0o5Zj0x/ozk7/J0+t55FNg13VZfnVqJlFxWsE2VpF8B20bEWne8KbmttQoYERGPN8K5GvV4WZIUJN1/dXU/mVnG3ILaxEjaTdIeabfPIOBE4O6C7YdK2lLJ7cjV112e3oDzNerxzMxKtdl887sZ6UjSrbc98B7wW5Kuqmr7kHRJtAFeBr5RrAtpHTT28czMSuIuPjMzyyV38ZmZWS5l1sXXpUuXKC8vz+r0ZmaWE9OnT18YEV1rr88sQZWXl1NZWdlwQTMz26xJKjpqjLv4zMwsl5ygzMwsl5ygzMwslxq8BiXpT8ARwHsR8eUi20Uy1Mq/kQwsOSrWnvvFzAyAFStWUFVVxbJlyxoubJuVsrIyunfvTuvWrUsqX8pNEuNJxgOra1y0w0nGndqFZPyvq1i3ccDMrBmpqqqiY8eOlJeXk/x/a81BRLBo0SKqqqro2bNnSfs0mKAiYoqk8nqKHAncmM4V9HQ6LM526YCRdfr0nVeY+Yv9SwqyqS3d5ZvsfcxPsg7DrFlYtmyZk1MzJIltttmGBQsWlLxPY1yD2oE1R3euStetRdLJkiolVeZlBIsey/9Fh3/e3XBBM2s0Tk7N07r+3Rvje1DFzlg0+0TEOGAcQEVFRfQ57/8a4fQbJi+tODMzW1NjtKCqWHNm0u4UmZnUzCwvTjjhBLp168aXv7zWfV91euKJJ5DEddd9Pm3Y888/jyR+85vf1Lvv1VdfzY031j+92YwZM5g8eXKDcVRWVnL66aeXFnQDxo8fz2mnndYox2oKjZGgJgHHp9M/DAYWN3T9ycwsS6NGjeLBBx9c5/369u3Lbbd9PgfjxIkT6devX4P7nXLKKRx//PH1lik1QVVUVHDFFVc0HOxmoMEEJelWkhkhd5VUJelESadIOiUtMhl4jWRmyWtIZsI0M8utAw44gK233nqd99txxx1ZtmwZ8+fPJyJ48MEHOfzww2u2X3PNNey1117069ePo446ik8++QSAsWPH1rSyhgwZws9+9jMGDRpEr169ePLJJ1m+fDljxozhtttuo3///tx22208++yz7LvvvgwYMIB9992X2bNnA0lL7ogjjqg57gknnMCQIUPYeeed10hcN998M4MGDaJ///58//vfZ9WqZMb566+/nl69enHggQcyderU9avAjaSUu/iOa2B7AP/VaBGZWbNx4b0zefntjxr1mL2378QFX++zzvtddtllTJgwYa31BxxwwBof/EcffTR33HEHAwYMYODAgbRt27Zm27e+9S1OOukkAEaPHs11113HD3/4w7WOuXLlSp599lkmT57MhRdeyKOPPspFF11EZWUlf/zjHwH46KOPmDJlCq1ateLRRx/lvPPO489//vNax3rllVd4/PHHWbJkCbvuuiunnnoqc+bM4bbbbmPq1Km0bt2aH/zgB0yYMIFDDjmECy64gOnTp9O5c2eGDh3KgAED1rmuNhZPWGhmBpx11lmcddZZDZb79re/zXe+8x1eeeUVjjvuOP7+97/XbHvppZcYPXo0H374IUuXLuXQQw8teoxvfetbAOy5557MnTu3aJnFixczcuRI/vnPfyKJFStWFC33ta99jbZt29K2bVu6devG/Pnzeeyxx5g+fTp77bUXAJ9++indunXjmWeeYciQIXTtmgwc/p3vfIdXX321wd85K05QZpaZ9WnpNJVSW1DbbrstrVu35pFHHuHyyy9fI0GNGjWKe+65h379+jF+/HieeOKJoueqbnW1bNmSlStXFi3z3//93wwdOpS7776buXPnMmTIkHqPVXi8iGDkyJH88pe/XKPsPffcs0nd4u8EBaxcFUya8VbWYdTo2K41Q3ftlnUYZs1KqS0ogIsuuoj33nuPli1brrF+yZIlbLfddqxYsYIJEyawww5FvxJaVMeOHVmyZEnN8uLFi2v2Hz9+fMnHATj44IM58sgjOfPMM+nWrRvvv/8+S5YsYe+99+aMM85g0aJFdOrUiTvuuKOkmzyy4gSV6tqxLOsQaixY4jHKzJrScccdxxNPPMHChQvp3r07F154ISeeeGLJ+++7775F11988cXsvffe7LTTTvTt23eNhNOQoUOHcumll9K/f3/OPfdczj77bEaOHMnvfvc7DjrooJKPA9C7d29+/vOfM2zYMFavXk3r1q258sorGTx4MGPHjmWfffZhu+22Y+DAgTU3T+SRshrRoaKiIvIwYeHMX+zPylXBJ9+dlHUoNRYsWcbw/qX/52W2KZk1axa777571mFYRor9/SVNj4iK2mU93YaZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZNSvz5s1j6NCh7L777vTp04fLL7+85H2rp9y49957a9YdccQRdY4YUe0Pf/hDzcCxja2u72Stz/m/973v8fLLLzdGWJSXl7Nw4cINOoYTlJk1K61ateK3v/0ts2bN4umnn+bKK69cpw/l7t27c8kll6zTOZsyQRUOtbSh57/22mvp3bt3Y4TVKJygzKxZqR5BAZLhhXbffXfeeqv0oc769etH586deeSRR9ba9thjjzFgwAD69u3LCSecwGeffcYVV1zB22+/zdChQxk6dOha+5SXl3Peeeexzz77UFFRwXPPPcehhx7KF7/4Ra6++moAli5dysEHH8zAgQPp27cvf/nLX2r279ChA5C07oYMGcLRRx/NbrvtxogRI4iIouc/9dRTqaiooE+fPlxwwQU1xxoyZAjVAyh06NCB888/n379+jF48GDmz58PwIIFCzjqqKPYa6+92GuvvWqm7Fi0aBHDhg1jwIABfP/736cxBoHwUEdmlp0HzoF3/9G4x9y2Lxx+aUlF586dy/PPP8/ee+8NlD5g7OjRoxk9ejSHHHJIzbply5YxatQoHnvsMXr16sXxxx/PVVddxY9+9CN+97vf8fjjj9OlS5eicfTo0YOnnnqKM888k1GjRjF16lSWLVtGnz59OOWUUygrK+Puu++mU6dOLFy4kMGDBzN8+PC1Bn59/vnnmTlzJttvvz377bcfU6dO5fTTT1/r/Jdccglbb701q1at4uCDD+bFF19kjz32WONYH3/8MYMHD+aSSy7h7LPP5pprrmH06NGcccYZnHnmmey///68+eabHHroocyaNYsLL7yQ/fffnzFjxnD//fczbty4kv4G9XGCMrNmaenSpRx11FH84Q9/oFOnTkDpA8Z+5StfAeDJJ5+sWTd79mx69uxJr169ABg5ciRXXnklP/rRjxo83vDhw4Fkxt6lS5fSsWNHOnbsSFlZGR9++CHt27fnvPPOY8qUKbRo0YK33nqL+fPns+22265xnEGDBtG9e3cA+vfvz9y5c9l///3XOt/tt9/OuHHjWLlyJe+88w4vv/zyWgmqTZs2NRMj7rnnnjUtxkcffXSNLtGPPvqIJUuWMGXKFO666y4gmQJkq622avD3bogTlJllp8SWTmNbsWIFRx11FCNGjKiZmwlKb0EBnH/++VxyySW0apV8jG5Il1b1lBktWrRYY/qMFi1asHLlSiZMmMCCBQuYPn06rVu3pry8nGXL1h5UutjUG7W9/vrr/OY3v2HatGlstdVWjBo1quixWrduXdNCKzzW6tWreeqpp2jXrt1a+zT2VB6+BgWsWLU66xDMbCOJCE488UR23313fvzjH6+x7ayzzmLGjBlrPWonJ4Bhw4bxwQcf8MILLwCw2267MXfuXObMmQPATTfdxIEHHgisPZXGulq8eDHdunWjdevWPP7447zxxhvrtH/h+T/66CPat29P586dmT9/Pg888MA6HWvYsGE1s/4CzJgxA0iSeHVyf+CBB/jggw/W6bjFOEEBK1dnM6K7mW18U6dO5aabbuKvf/0r/fv3p3///kyePHm9jnX++edTVVUFQFlZGddffz3HHHMMffv2pUWLFpxyyikAnHzyyRx++OFFb5IoxYgRI6isrKSiooIJEyaw2267rdP+hefv168fAwYMoE+fPpxwwgnst99+63SsK664gsrKSvbYYw969+5dcyPHBRdcwJQpUxg4cCAPP/wwO+644zodtxhPt/GL/VmybCUx8v6sQ6nh6TZsc+bpNpo3T7dhZmabPCcoMzPLpZLu4pN0GHA50BK4NiIurbV9R+AGYMu0zDkRsX6duhm56L6ZWYdQY48dOruLzzZrEdHod3xZ/q3rJaUGE5SklsCVwCFAFTBN0qSIKBwbZDRwe0RcJak3MBkoX6dINoLfP/Iqlz/2zzXWTWyT3Do5652177Dp0qENXTu2XWt9U3pj0SesWOm7Cm3zVVZWxqJFi9hmm22cpJqRiGDRokWUlZWVvE8pLahBwJyIeA1A0kTgSKAwQQXQKX3eGXi75Ag2ojMP6cWZh/Rac+X1/8PTry/i1pMGZxNULRfdN9MJyjZr3bt3p6qqigULFmQdim1kZWVlNV8kLkUpCWoHYF7BchWwd60yY4GHJf0QaA98teQIzKxZad26NT179sw6DNsElJKgirXBa3ckHgeMj4jfStoHuEnSlyNijaaApJOBk4FGuUe+MW3z+r0NF9oIWn/aGcL3rpiZlZKgqoAeBcvdWbsL70TgMICIeEpSGdAFeK+wUESMA8ZB8j2o9Yy5Saxo1zXrEACIlqAVn2UdhplZ5kr5V30asIuknpLaAMcCk2qVeRM4GEDS7kAZ4A5mMzNbbw0mqIhYCZwGPATMIrlbb6akiyQNT4v9BDhJ0gvArcCoyGqIivXQte2qrEMwM7NaSvoeVPqdpsm11o0peP4ysG4DOuVIt7YreTPrIMzMbA2+Gm9mZrnk+aByKAgmzSh9Cuqm1LFda4bu2i3rMMysGXKCyqHWLVrQtWPp37ZuSguWrD2RmZnZxuAuPjMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzy6WSEpSkwyTNljRH0jl1lPm2pJclzZR0S+OG2bTaLP8g6xDMzKyWVg0VkNQSuBI4BKgCpkmaFBEvF5TZBTgX2C8iPpDUrakCbgptnaDMzHKnlBbUIGBORLwWEcuBicCRtcqcBFwZER8ARMR7jRummZk1N6UkqB2AeQXLVem6Qr2AXpKmSnpa0mHFDiTpZEmVkioXLFiwfhGbmVmzUEqCUpF1UWu5FbALMAQ4DrhW0pZr7RQxLiIqIqKia9eu6xqrmZk1Iw1egyJpMfUoWO4OvF2kzNMRsQJ4XdJskoQ1rVGi3AjKKy/OOgQAfvEpTGmxN3BA1qGYmWWqlAQ1DdhFUk/gLeBY4Lu1ytxD0nIaL6kLSZffa40ZaKN4/Jfwt0uLbmr/way11i0v68KKdhu3pddz9RsQq1m6Uc9qZpY/DSaoiFgp6TTgIaAl8KeImCnpIqAyIial24ZJehlYBZwVEYuaMvD1MvTc5FHb2M7MPCQfd8avfuxiiNVZh2FmlrlSWlBExGRgcq11YwqeB/Dj9GFmZrbBSkpQtvFt8/q9WYeQWNmWpPfWzGzjcoLKI7XY6Ne+6rRwHpNmvJV1FGa2GWvRtkOnYuudoHJo+WooyzqI1Dbty6BjXqIxs81SixYti67e2HHk0dvbHZJ1CGtY4XskzMycoADe3v7QrEMwM7NanKDMzCyXnKDMzCyXfJNETp3z96wjSBzUtYzBPbOOwsyaIyeoDE2YDbe8uua6iW2Sn/8oMg5Ht3bwhS2aPq5qry0GrWrL4I13SjOzGk5QGRqxa/IoVF6ZJKf7v55NTIXO+TvJwFVmZhnwNSgzM8slJygzM8slJygzM8slJygzM8slJ6gc6tYu6wjMzLLnBJVDG/NWcjOzvHKCMjOzXHKCMjOzXPIXdYF2bVry4SfLsw4DgJWrVrM6so7CzCx7TlBAn+07QYcuWYeRmNOaxZ+uyDoKM7PMuYvPzMxyyQnKzMxyqaQEJekwSbMlzZF0Tj3ljpYUkioaL0QzM2uOGkxQkloCVwKHA72B4yT1LlKuI3A68ExjB2lmZs1PKS2oQcCciHgtIpYDE4Eji5S7GPg1sKwR4zMzs2aqlAS1AzCvYLkqXVdD0gCgR0Tc14ixmZlZM1ZKglKRdTXf1JHUAvg98JMGDySdLKlSUuWCBQtKj9LMzJqdUhJUFdCjYLk78HbBckfgy8ATkuYCg4FJxW6UiIhxEVERERVdu3Zd/6jNzGyzV0qCmgbsIqmnpDbAscCk6o0RsTgiukREeUSUA08DwyOiskkiNjOzZqHBBBURK4HTgIeAWcDtETFT0kWShjd1gGZm1jyVNNRRREwGJtdaN6aOskM2PCwzM2vuPJKEmZnlkhOUmZnlkhOUmZnlkhOUmZnlkhOUmZnlkhOUmZnlkhOUmZnlkhOUmZnlkhOUmZnlkhOUmZnlkhOUmZnlUklj8dnG1f7jNymvvDjrMPjFp3Dv6sHAAVmHYmabsRZlHbYqtt4JKm92HsLHn67IxR+m5+o3OGh1weyUZmZNoEXbLZygNgm9DuOlFhVsuUWbrCNh9WMXw+qsozCz5srXoMzMLJfcggIo6wxL52cdRY0tPn4Httg56zDMzDLlBAWwyyFZR7CGVXOvyTqENVx038ysQzCzZsgJygCYMBtueXXNdRPTy2Cz3lmyVvkuHdrQtWPbjRCZmW0uFiz5jIVLl5dc3gnKABixa/IoVF4J/1gEt540OJugzKxZGHJD8fW+ScLMzHLJCcrMzHLJCcrMzHLJCcrMzHKppAQl6TBJsyXNkXROke0/lvSypBclPSZpp8YP1bLwhbJVWYdgZpu51Z998kGx9Q0mKEktgSuBw4HewHGSetcq9jxQERF7AHcCv96wcC0vvlDmsY7MrGmtXrZ0/RIUMAiYExGvRcRyYCJwZGGBiHg8Ij5JF58Gum9IsGZmZqUkqB2AeQXLVem6upwIPFBsg6STJVVKqlywYEHpUZqZWbNTyhd1VWRd0RkYJP07UAEcWGx7RIwDxgFUVFR4Foc6tGvTkg8/Kf3b1k1l5arVsNpdfGaWjVISVBXQo2C5O/B27UKSvgqcDxwYEZ81TnjNU5/tO0GHLlmHAXNas+TjFVlHYWbNVCldfNOAXST1lNQGOBaYVFhA0gDgf4HhEfFe44dpZmbNTYMJKiJWAqcBDwGzgNsjYqakiyQNT4tdBnQA7pA0Q9KkOg5nZmZWkpIGi42IycDkWuvGFDz/aiPHZWZmzZxHkjAzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1wqaboN28jKOsPS+VlHAatWEGqZdRRm1kw5QeXRLodkHUGi8npi+YKsozCzZspdfGZmlktOUGZmlktOUGZmlku+BmX1atFCLFiyLOswzGxztnr1qmKrnaCsXh3btmJ4/x2yDsPMNmOrP1v6UbH17uIzM7NcKilBSTpM0mxJcySdU2R7W0m3pdufkVTe2IGamVnz0mCCktQSuBI4HOgNHCepd61iJwIfRMSXgN8Dv2rsQM3MrHkppQU1CJgTEa9FxHJgInBkrTJHAjekz+8EDpakxuSzPXMAAA8HSURBVAvTzMyam1IS1A7AvILlqnRd0TIRsRJYDGxT+0CSTpZUKalywQKPUJB72/aFrXfOOgoza6ZKuYuvWEso1qMMETEOGAdQUVGx1nbLmcMvzToCM2vGSmlBVQE9Cpa7A2/XVUZSK6Az8H5jBGhmZs1TKQlqGrCLpJ6S2gDHApNqlZkEjEyfHw38NSLcQjIzs/XWYBdfRKyUdBrwENAS+FNEzJR0EVAZEZOA64CbJM0haTkd25RBm5nZ5q+kkSQiYjIwuda6MQXPlwHHNG5oZmbWnHkkCTMzyyUnKDMzyyVldS+DpCXA7ExOvn66AAuzDmIdON6m5XibluNtWnmLd6eI6Fp7ZZajmc+OiIoMz79OJFU63qbjeJuW421ajrdpuIvPzMxyyQnKzMxyKcsENS7Dc68Px9u0HG/TcrxNy/E2gcxukjAzM6uPu/jMzCyXnKDMzCyXMklQDU0hnzVJPSQ9LmmWpJmSzkjXby3pEUn/TH9ulXWs1SS1lPS8pPvS5Z6SnkljvS0d6Dc3JG0p6U5Jr6T1vE/O6/fM9LXwkqRbJZXlqY4l/UnSe5JeKlhXtD6VuCJ9/70oaWBO4r0sfT28KOluSVsWbDs3jXe2pEPzEG/Btp9KCkld0uVc1m+6/odpHc6U9OuC9ZnWb50iYqM+SAac/RewM9AGeAHovbHjaCDG7YCB6fOOwKsk093/GjgnXX8O8KusYy2I+cfALcB96fLtwLHp86uBU7OOsVa8NwDfS5+3AbbMa/2STMj5OtCuoG5H5amOgQOAgcBLBeuK1ifwb8ADJPO4DQaeyUm8w4BW6fNfFcTbO/2caAv0TD8/WmYdb7q+B8lA2m8AXXJev0OBR4G26XK3vNRvXY8sWlClTCGfqYh4JyKeS58vAWaRfEgVTm1/A/CNbCJck6TuwNeAa9NlAQcBd6ZFchMrgKROJG+g6wAiYnlEfEhO6zfVCmiXzne2BfAOOarjiJjC2nOw1VWfRwI3RuJpYEtJ222cSBPF4o2IhyOZkRvgaZK55yCJd2JEfBYRrwNzSD5HNpo66hfg98DZrDlBay7rFzgVuDQiPkvLvJeuz7x+65JFgiplCvnckFQODACeAb4QEe9AksSAbtlFtoY/kLxJVqfL2wAfFrzZ81bHOwMLgOvTbslrJbUnp/UbEW8BvwHeJElMi4Hp5LuOoe763BTegyeQtEIgp/FKGg68FREv1NqUy3iBXsBX0m7pv0naK12f13gzSVAlTQ+fB5I6AH8GfhQRH2UdTzGSjgDei4jphauLFM1THbci6X64KiIGAB+TdEHlUnrt5kiS7o/tgfbA4UWK5qmO65Pr14ek84GVwITqVUWKZRqvpC2A84ExxTYXWZeH+m0FbEXS7XgWcHva25LXeDNJUKVMIZ85Sa1JktOEiLgrXT2/uqme/nyvrv03ov2A4ZLmknSXHkTSotoy7Y6C/NVxFVAVEc+ky3eSJKw81i/AV4HXI2JBRKwA7gL2Jd91DHXXZ27fg5JGAkcAIyK9QEI+4/0iyT8sL6Tvve7Ac5K2JZ/xQhLXXWnX47MkPS5dyG+8mSSoUqaQz1T6X8V1wKyI+F3BpsKp7UcCf9nYsdUWEedGRPeIKCepy79GxAjgceDotFguYq0WEe8C8yTtmq46GHiZHNZv6k1gsKQt0tdGdby5reNUXfU5CTg+vdtsMLC4uiswS5IOA34GDI+ITwo2TQKOldRWUk9gF+DZLGKsFhH/iIhuEVGevveqSG6sepec1i9wD8k/sEjqRXJz0kJyWL81srgzg+Qul1dJ7hY5P8u7ROqIb3+SJu6LwIz08W8k13YeA/6Z/tw661hrxT2Ez+/i25nkRTYHuIP0zp28PID+QGVax/eQdD3ktn6BC4FXgJeAm0jueMpNHQO3klwfW0HyYXliXfVJ0qVzZfr++wdQkZN455BcC6l+z11dUP78NN7ZwOF5iLfW9rl8fhdfXuu3DXBz+hp+DjgoL/Vb18NDHZmZWS55JAkzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJyjbLCmZzuMHBcvbS7qzvn3W4dhjJf00fX6RpK820nG3UzpdSlOQtHQdyj6qHE13Ys2TE5RtrrYEahJURLwdEUfXU369RMSYiHi0kQ73Y+CaRjrWhrqJgvozy4ITlG2uLgW+KGlGOhFeefXkbZJGSbpH0r2SXpd0mqQfpyOrPy1p67TcFyU9KGm6pCcl7Vb7JJLGSzo6fT5X0oWSnpP0j+ryktqnE8hNS89R1/QyRwEPpvtMlrRH+vx5SWPS5xdL+l76/Kz0mC9KurAgpn+X9Gz6u/+vpJa1Yu4i6SlJX0tbbVPSsi9J+kpabBJw3HrWvVmjcIKyzdU5wL8ion9EnFVk+5eB75LMe3MJ8EkkI6s/BRyflhkH/DAi9gR+CvxPCeddGBEDgavSfSAZRuavEbEXyaRxl6XTi9RIx0D7INK5eoApJFMjdCIZ2Xu/dP3+wJOShpGMmTaIZNioPSUdIGl34DvAfhHRH1gFjCg4zxeA+4ExEXF/WgcPpWX7kQwxRER8ALSVtE0Jv7NZk2jVcBGzzdLjkUxGuUTSYuDedP0/gD3SqVb2Be5IxocFkvH3GlI98v104Fvp82EkI85XJ6wyYEeSiTCrbUcyR1a1J4HTSWbyvR84JJ3ioTwiZks6KT3u82n5DiQJaw9gT2BaGnc7Ph/FvDXJmHz/FRF/S9dNA/6Ujt5/T0TMKIjhPZLpRRaV8HubNTonKGuuPit4vrpgeTXJ+6IFyYSE/dfzuKv4/P0l4KiImF3Pfp+SJK5q04AK4DXgEZJpEU4iSXzVx/xlRPxv4UEk/RC4ISLOLXKOlen+hwJ/g2TmVUkHkMzIfJOkyyLixrR8WRqXWSbcxWebqyVAx/XdOZIJKl+XdAwkU7BI6reeh3sI+GE6VQeSBhQp8ypQXnD+5SQje3+bZPrzJ0m6DJ8sOOYJaUsPSTtI6kbSQjo6fY6krSXtVH1Ykplqd5N0Trp9J5IJL68hmWJmYPXvC2xLMkq3WSacoGyzFBGLgKnphf/L1vMwI4ATJb0AzCSZVXd9XEzSvfZieqPGxUXi/Rj4l6QvFax+EpgfydxIT5JMJPdkWv5h4BbgKUn/IJn0sWNEvAyMBh6W9CJJ62u7gvOsIpk3bGh6G/4QYIak50lu0rg8Lbon8HR8PqW92Ubn6TbMckLSN4E9I2J0DmK5HJgUEY9lHYs1XxvcgpL0TUlReAtuekvvd9fhGH9vYHvJXzA021RFxN3kp0vtJScny9oGt6Ak3U7ShfBYRIxN1w0BfhoRRzSwb8u0y6GhcyyNiA4bFKiZmW1SNqgFlV6g3Y9kOuFjCzZdSvIdjhmSzqy1zxBJj0u6heSW3poWUj1fGqzet+YLhhsSt5mZ5d+G3mb+DeDBiHhV0vuSBkbEcyRfkqyvBTUI+HJEvF5rffWXBi9Jv/2+RfWG9AuGk4DREfHIBsZtZmY5t6HXoI4DJqbPJ1L60CjPFklOkHz34z8ljQX6pl+khM+/YHi2k5OZWfOw3gkqHQLlIOBaSXOBs4DvVH/XowEfF1sZEVOAA4C3SL40WD3kTOEXDM3MrBnYkBbU0cCNEbFTRJRHRA+SYVn2Zz2/JFnXlwYp8gVDMzPbvG1IgjoOuLvWuj+TXEd6EVgp6YXaN0k0YAjFvzRY7AuGZma2GfMXdc3MLJc81JGZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeVSoyUoSX+S9F46IVv1umMkzZS0WlJFY52rqUk6M437JUm3SipreK9sSWop6XlJ92UdS0OKvVbyzPE2LcfbtDa1eAs1ZgtqPHBYrXUvAd8CpjTieZqUpB2A04GKiPgy0JI1R2rPqzOAWVkHUaLxrP1aybPxON6mNB7H25TGs2nFW6PRElQ6jt77tdbNiojZjXWOjagV0E5SK5IR1d/OOJ56SeoOfA24NutYSlHstZJnjrdpOd6mtanFW8jXoGqJiLeA3wBvAu8AiyPi4WyjatAfgLOB1VkHYmbWWJygapG0FXAk0BPYHmgv6d+zjapuko4gGWB3etaxmJk1JieotX0VeD0iFkTECuAuYN+MY6rPfsDwdMqTicBBkm7ONiQzsw3nBLW2N4HBkrZI57Y6mBzffBAR50ZE94goJ7mZ468RkdsWn5lZqRrzNvNbgaeAXSVVSTpR0jclVQH7APdLeqixztdUIuIZ4E7gOeAfJHU0LtOgNjPFXitZx1Qfx9u0HG/T2tTiLeTpNszMLJfcxWdmZrnkBGVmZrnkBGVmZrnU1GPxXSbpFUkvSrpb0paNdb6mJGlLSXemsc+StE/WMdVHUpmkZyW9kI4heGHWMTVE0mGSZkuaI+mcrONpiONtWo63aW1q8daIiEZ5AAcAA4GXCtYNA1qlz38F/KqxzteUD+AG4Hvp8zbAllnH1EC8Ajqkz1sDzwCDs46rnnhbAv8Cdk7r9wWgd9ZxOV7H63jz9WjqsfgejoiV6eLTQPfGOl9TkdSJJNleBxARyyPiw2yjql8klqaLrdNHnm/PHATMiYjXImI5yReMj8w4pvo43qbleJvWphZvjY15DeoE4IGNeL71tTOwALg+nb7iWkntsw6qIel0GzOA94BHIvk+V17tAMwrWK5K1+WV421ajrdpbWrx1tgoCUrS+cBKYMLGON8GakXSVXlVRAwAPgZy32cbEasioj9JK3WQpC9nHVM9VGRdnlt8jrdpOd6mtanFW6PJE5SkkcARwIhIO0RzrgqoKmiB3EmSsDYJaXfkE+R7/pcqoEfBcnfyPaWJ421ajrdpbWrx1mjSBCXpMOBnwPCI+KQpz9VYIuJdYJ6kXdNVBwMvZxhSgyR1rb5DUlI7kgFvX8k2qnpNA3aR1FNSG5IxBCdlHFN9HG/TcrxNa1OLt0arxjpQOt7TEKBLOv7eBcC5QFvgkWTcVZ6OiFMa65xN6IfAhPSP+RrwnxnH05DtgBsktST5p+P2iMjt1O8RsVLSacBDJHcY/SkiZmYcVp0cb9NyvE1rU4u3kMfiMzOzXPJIEmZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlkv/HwaiPnd2zQ5aAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "stride = 16\n",
    "agg_d1, agg_n1 = aggregate(d1, n1, stride)\n",
    "agg_d2, agg_n2 = aggregate(d2, n2, stride)\n",
    "agg_d1, agg_n1, agg_d2, agg_n2 = [list(map(int, mpc.run(mpc.output(_)))) for _ in (agg_d1, agg_n1, agg_d2, agg_n2)]\n",
    "T1_, E1_ = events_from_table(agg_d1, agg_n1)\n",
    "T2_, E2_ = events_from_table(agg_d2, agg_n2)\n",
    "T1_, T2_ = [t * stride for t in T1_], [t * stride for t in T2_]\n",
    "fit_plot(T1_, T2_, E1_, E2_, 'Aggregated Kaplan-Meier survival curves', 'weeks', '1=Maintained', '2=Not maintained')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Picking `stride = 16` achieves a reasonable balance between privacy and utility. To enhance both privacy and utility at the same time, one may look for differentially private randomization techniques, adding a suitable type of noise to the Kaplan-Meier curves. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Secure Logrank Tests\n",
    "\n",
    "The function `logrank_test()` performs a secure logrank test on a secret-shared dataset, similar to function `lifelines.statistics.logrank_test()` used above for a dataset in the clear. The input parameter `secfxp` specifies the secure type to be used for fixed-point arithmetic, and the output is an instance of `lifelines.statistics.StatisticalResult`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.06533935619386949\n"
     ]
    }
   ],
   "source": [
    "print((await logrank_test(secfxp, d1, d2, n1, n2)).p_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Relying solely on p-values is in general not a good idea, and this is especially true when handling otherwise (mostly) hidden data. Together with the aggregated curves, however, the p-value may lead to a useful conclusion for a study. \n",
    "\n",
    "The function `logrank_test()` uses one secure fixed-point division per time moment in $1..maxT$. Even though these divisions can all be done in parallel, the total effort is significant when $maxT$ is large. However, \"most of the time\" there is actually no event happening and no divisions need to be performed for these time moments. E.g., in the survival tables for the aml dataset above, there are only 7 time moments with nonzero `d1` entries on the entire timeline $1..161$, and only 9 time moments with nonzero `d2` entries.\n",
    "\n",
    "Therefore, it may be advantageous to first extract the nonzero rows of the survival tables, and then limit the computation of the logrank test to those rows. The extraction needs to be done obliviously, not leaking any information about (the location of) the nonzero entries of the survival tables. To prevent this oblivious extraction step from becoming a bottleneck, however, we will actually exploit the fact that the aggregate curves are revealed anyway. We may simply use `agg_d1` and `agg_d2` to bound the number of events per stride, and extract the nonzero rows obliviously and efficiently for each stride. \n",
    "This is basically what the function `agg_logrank_test()` does:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.06533931091264207\n"
     ]
    }
   ],
   "source": [
    "print((await agg_logrank_test(secfxp, d1, d2, n1, n2, agg_d1, agg_d2, stride)).p_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even for a small dataset like aml, the speedup is already noticeable. For larger datasets, the speedup gets really substantial, as can be noticed for some of the other datasets included with the [kmsurival.py](kmsurvival.py) demo. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "We end with two complete runs of the demo on the aml dataset, showing the Chi2 test statistic and p-value for each logrank test. \n",
    "\n",
    "The help message included with the demo shows the command line options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Showing help message for kmsurvival.py, if available:\n",
      "\n",
      "usage: kmsurvival.py [-h] [-i I] [-s S] [-a A] [--collapse] [--print-tables]\n",
      "                     [--plot-curves]\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help          show this help message and exit\n",
      "  -i I, --dataset I   dataset 0=btrial(default) 1=waltons 2=aml 3=lung 4=dd\n",
      "                      5=stanford_heart_transplants 6=kidney_transplant\n",
      "  -s S, --stride S    interval length for aggregated events\n",
      "  -a A, --accuracy A  number of fractional bits\n",
      "  --collapse          days->weeks->month->years\n",
      "  --print-tables      print survival tables\n",
      "  --plot-curves       plot survival curves\n"
     ]
    }
   ],
   "source": [
    "!python kmsurvival.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Complete Run:  5 logrank tests + survival curves\n",
    "To show the plots of the survival curves the `main()` function of the demo is called directly from a notebook cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using secure fixed-point numbers: SecFxp64:32\n",
      "Dataset: aml, with 1-party split, time 1 to 161 (stride 16) weeks\n",
      "Chi2=3.396389, p=0.065339 for all events in the clear\n",
      "Chi2=3.396389, p=0.065339 for own events in the clear\n",
      "Chi2=2.685357, p=0.101275 for aggregated events in the clear\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEOCAYAAADc94MzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwU1b3//9ebXQRBATdAQa8bOziumIgSF6JiFpKI3iiSaGKixptEQ9QkLr/ca6I3uSZfs7gSjWIMiUoSjEuCcY+4oIKIQZYIEgQiCLII+Pn9UTVjM/TM9Mz0TFfPvJ+PRz+mu+pU1adrevoz59SpcxQRmJmZZU2bUgdgZmaWjxOUmZllkhOUmZllkhOUmZllkhOUmZllkhOUmZllkhOUlQ1JiyR9LANxnCHpoSLsZ4KkJxqx/RxJoxq47WRJ/1/6/COS5hW43QOSzqphXT9JIaldQ2KqKT5rvZygMij9It4gaZ2k5ZJuk9SlgfsaJWlJI+P5L0n/krRG0q2SOjZmf01BUh9Jv5O0Mo3zFUkTmuJYEXFnRBzfFPuuZxwDI+LRIuzn8Yg4oMCyYyLiV409ZjnKyj9IrYkTVHadEhFdgBHAIcDl9d1Bkf6TPQGYBIwG+gH7AFc2dr9N4A7gTWBvoAdwJrC8ITsqxnkzs8Zzgsq4iFgKPAAMApB0tqS5ktZKWiDpS5VlK2tLkr4l6V/AlHTbPdPa2DpJe0paL6lHznYHS1ohqX2eEM4CbomIORHxDnA1MKHQ+CX9Nqf29ZikgTnrJkv6WdpstE7Sk5J2l/R/kt6R9Jqk4QUe6hBgckS8FxFbIuLFiHgg97xUi6vqv2FJV0iaKunXkt4FLk1rsLvklB+e1s7a5zbNSfqFpOuq7ft+SV9Pn0+S9Eb6+3pV0icLPXd1yfMe7pF0e3qsOZIqqsX/QrruN0CnnHVV5yeNd2q141wv6Sfp80clfTF93lbSdel5WQCcVFN8OTH+Oud1jZ+NAt77xPTv4B1JD0raO11e1+9jz7SmvULSQkkXVosv7zmUdAewF/CH9LN6iaRO6WdmlaTVkmZK2q3Q92B1c4LKOEl9gY8DL6aL3gZOBnYCzgZ+LGlEzia7A7uQ1CTOBMYAb0VEl/TxFvAo8Nmcbf4TuDsiNucJYSDwUs7rl4DdKhOcpD9KmlTLW3gA2A/YFXgBuLPa+s+S1A57ApuAp9NyPYGpwI9q2XeuZ4AbJJ0maa8Ct8l1anq87sC1aRyfzll/OjA1zzm6C/icJAFI2hk4Hrg7Xf8G8BGgG0nN89eS9mhAfIUYmx63OzAN+H9pTB2A+0hqmbsAv2Xb95ZrCvBxSTul27Yl+R3dlafsOSSfxeFABTCunvHW9dnIS9IngEuBTwG9gMfTuKGW34ekNsAfSD7DvUlaBS5S0kpQKe85jIjPA/8kbdmIiB+S/PPWDehLUmv/MrChnufAauEElV33SVoNPAH8DfhvgIj4U0S8EYm/AQ+RfAFW+gD4XkRsioia/lh+RZKUKr+AxpN8eeXTBViT87ryedc0npMj4pqa3kRE3BoRayNiE3AFMFRSt5wi90bE8xGxEbgX2BgRt0fEVuA3JF9+hfgMyRfVd4CFkmZJOqTAbQGejoj7IuKD9LzdRXJeSL/sTiP/l/TjQPDh72Bcuq+3ACLitxHxVrrf3wD/AA6tR1z18URETE/P3R3A0HT54UB74P8iYnNETAVm5ttBRCwmSRafSBcdC6yPiGfyFP9sus83I+LfwP/UJ9gCPhs1+RLwPxExNyK2kPxtDEtrUbX9Pg4BekXEVRHxfkQsAG4i+d1Wqukc5rOZJDH9R0RsTT/H7xb6/q1uTlDZ9YmI6B4Re0fEVyqTjaQxkp6R9O80gX2cpLZRaUX6ZV+b+4EBkvYBjgPWRMSzNZRdR1Jbq1T5fG1dbyBtArombeJ6F1iUrsqNN/c60YY8rwvqHBIR70TEpIgYCOwGzCJJ8ipke5LrV7mmAkdI2hP4KMmX3uN5jhsk/3GPTxedTk5NQNKZabJcnf6+BrHt+88rbaqqbJa9tMD38K+c5+uBTkqup+0JLI1tR4ZeXMt+qpJz+n7yJWbS/eaet9r2uY0CPxs12Ru4Puec/hsQ0LuO38feJM3dq3O2vZTk81KppnOYzx3AgyS1s7ck/VD5m8mtgZygyoiS3nO/A64DdouI7sB0kj/OStWHp99uuPo0gd0DnAF8npprTwBz2Pa/yKHA8ohYVUDIp5M0nX2MpCmkX+VbKWDbBouIlSTnaE+SJq33gM6V69NaY6/qm1Xbx2qS2ulnSd7HlGpf8LmmAOPS/+API/kdkb6+CTgf6JH+vmZTwPuPiC/nNMv+d13l67AM6F0tWdfWDPpbYJSkPsAnqTlBLSNp3qppn9ucd5Lm50qN+Wy8CXwp/Qeu8rFDRDyVrs/7+0i3W1htu64R8fECjgnbf0Y2R8SVETEAOJKkufPMAvdlBXCCKi8dgI7ACmCLpDEk7eu1WQ70yNN0cjtJZ4exwK+rb1St3BckDUjb8y8HJhcYb1eS60qrSL6oGvtFWyNJP5A0SFI7SV2B84D5aSJ9neQ/4ZPS/3AvJzmPdbmL5Avn09T8JU1EvEjyO7kZeDBNbgA7knyprUhjPJu0s0szexrYAlyYnp9PUUszY0SsILlOeRvJF/rcGorek+6zT/rZqH4tchZwmpKOJdWvUTXms/EL4NuVnSokdZP0mZz4a/p9PAu8q6QT0Q5pLW5QPZqCl5P0YiU97jGSBqf/8LxL0uS3tR7vw+rgBFVGImItcCHJF8M7JP+FTqtjm9dI/qNckDZr7Jkuf5LketULEbGolu3/DPwQmEHShLMY+F7leiU98Gpqgro9Lb8UeJWkI0NT6UxyDWs1sICkOWcsQESsAb5C8oW1lOQ/+0LuDZtGchF/eUS8VEfZKSS1gapEFhGvAv9LkiCWA4OBJwt+R0USEe+TdCiYQPK5+Rzw+zo2u4tq7yePm0iauF4iuW5VfZ/fAfZNj3lltX01+LMREfcCPyBpWnuXpFY6plqxfL+PrcApwDBgIbCS5DNRyHUvSK6xXZ7+HX2TpEY4lSQ5zSW5VlzbP3tWT6q51cJaOkl/Be6KiJtLHYuZWXVOUK1U2qzxMNA3rZmZmWWKm/haIUm/Ah4BLnJyMrOscg3KzMwyyTUoMzPLpJINitmzZ8/o169fqQ5vZmbN5Pnnn18ZEdXvPaxTyRJUv379eO6550p1eDMzayaSCh5lJJeb+MzMLJOcoMzMLJOcoMzMLJPqvAYl6VaSQRDfjojtxhFLB6C8nmRU7fXAhIh4odiBmlnLt3nzZpYsWcLGjXUNyG9Z1KlTJ/r06UP79sUZ1L2QThKTSSbtur2G9WNIxivbj2Tk4J+nP83M6mXJkiV07dqVfv36UfhMKZYFEcGqVatYsmQJ/fv3L8o+60xQEfGYpH61FDkVuD2diuAZSd0l7RERy2rb74ZlrzHnv48qONAndziGv3SufVT8U4f15vTDGjKZqpllwcaNG52cypQkevTowYoVK4q2z2Jcg+rNtpOWLUmXbUfSuZKek/RcfUaw6Ld5ASM3zKi1zKvL3uX+WUsL3qeZZZOTU/kq9u+uGPdB5Ysob/aJiBuBGwEqKipi4KVPFHaE205iIPCbs4+oscjnfvl0YfsyM7OyUIwa1BK2nVWzD/BWEfZbb6s3bGbarKV1PmbMe7sU4ZlZGZDE5z//+arXW7ZsoVevXpx88sm1bvfcc89x4YUX1lpm9erV/OxnPysojiOPPLKgcnVZtGgRgwaVYp7MxitGDWoacL6ku0k6R6yp6/pTU9m69QN6de1UZ7kVa91DyMzy23HHHZk9ezYbNmxghx124OGHH6Z377xXLbZRUVFBRUVFrWUqE9RXvvKVOvf31FNP1VmmpauzBiVpCsmMoAdIWiLpC5K+LOnLaZHpJDOYzieZYbPuM29mlmFjxozhT3/6EwBTpkxh/PjxVeueffZZjjzySIYPH86RRx7JvHnzAHj00UerallXXHEFEydOZNSoUeyzzz785Cc/AWDSpEm88cYbDBs2jIsvvph169YxevRoRowYweDBg7n//vurjtOlS5eq/Y4aNYpx48Zx4IEHcsYZZ1B5Df/555/n6KOP5uCDD+aEE05g2bJlVcuHDh3KEUccwQ033NDEZ6vpFNKLb3wd6wP4atEiMjMDrvzDHF59692i7nPAnjvxvVMG1lnutNNO46qrruLkk0/m5ZdfZuLEiTz++OMAHHjggTz22GO0a9eORx55hEsvvZTf/e532+3jtddeY8aMGaxdu5YDDjiA8847j2uuuYbZs2cza9YsIGk+vPfee9lpp51YuXIlhx9+OGPHjt2us8GLL77InDlz2HPPPRk5ciRPPvkkhx12GBdccAH3338/vXr14je/+Q2XXXYZt956K2effTY//elPOfroo7n44ouLcOZKo2SDxZqZZdWQIUNYtGgRU6ZM4eMf3/b2ljVr1nDWWWfxj3/8A0ls3rw57z5OOukkOnbsSMeOHdl1111Zvnz5dmUigksvvZTHHnuMNm3asHTpUpYvX87uu+++TblDDz2UPn36ADBs2DAWLVpE9+7dmT17NscddxwAW7duZY899mDNmjWsXr2ao48+GoDPf/7zPPDAA40+J6VQPglqw7/hlan513XqBnRp1nDMrGkVUtNpSmPHjuWb3/wmjz76KKtWrapa/p3vfIdjjjmGe++9l0WLFjFq1Ki823fs2LHqedu2bdmyZct2Ze68805WrFjB888/T/v27enXr1/eUTTy7SsiGDhwIE8/vW0P5tWrV7eYrvrlk6A+2Apddsu/bt1y6pOgNm7ZyrR63DPVdYf2HHPArgWXN7PyN3HiRLp168bgwYN59NFHq5avWbOmqtPE5MmT67XPrl27snbt2m32teuuu9K+fXtmzJjB4sWFz0pxwAEHsGLFCp5++mmOOOIINm/ezOuvv87AgQPp1q0bTzzxBEcddRR33nlnvWLMkvJJUEXUd+cd61Xevf7MWp8+ffrwta99bbvll1xyCWeddRY/+tGPOPbYY+u1zx49ejBy5EgGDRrEmDFj+Na3vsUpp5xCRUUFw4YN48ADDyx4Xx06dGDq1KlceOGFrFmzhi1btnDRRRcxcOBAbrvtNiZOnEjnzp054YQT6hVjlqg+IzoUU0VFRRQ8YeFtJ8F7K+Ck/82/ft1yPvdUb1at28TVnxhcvCBTK9ZuZOywuruZmlnjzJ07l4MOOqjUYVgj5PsdSno+Imrvg5+Hp9swM7NMcoIyM7NMcoIyM7NMalGdJJau3sBVf5xTZ7mR+/Zk9EE19Ag0M7NMaDE1qFOH9aZ39x3qLLd41XqefGNlM0RkZmaN0WJqUKcfthddOratc7DYQmpY1dXnvinfM2VmVhwtpgbVlPruvCO9unYq6LF2Q/5hT8ysPEjiG9/4RtXr6667jiuuuKLWbe677z5effXVJonni1/8Yp37LvT4v/jFL7j99tuLEteECROYOrWG0X2KxAnKzCxHx44d+f3vf8/KlYVfCmjKBHXzzTczYMCAohz/y1/+MmeeeWaxQmtyTlBmZjnatWvHueeey49//OPt1i1evJjRo0czZMgQRo8ezT//+U+eeuoppk2bxsUXX8ywYcN44403ttlmwoQJnHfeeRxzzDHss88+/O1vf2PixIkcdNBBTJgwoarceeedR0VFBQMHDuR73/te1fJRo0ZROahBly5duOyyyxg6dCiHH344y5cvz3v8m266iUMOOYShQ4fy6U9/mvXr1wPJNCDXXXdd1X6/9a1vceihh7L//vtXjda+detWLr74Yg455BCGDBnCL3/5SyAZ2Pb8889nwIABnHTSSbz9dtNP/NpirkGZWQvzwCT41yvF3efug2HMNXUW++pXv8qQIUO45JJLtll+/vnnc+aZZ3LWWWdx6623cuGFF3LfffcxduxYTj75ZMaNG5d3f++88w5//etfmTZtGqeccgpPPvkkN998M4cccgizZs1i2LBhfP/732eXXXZh69atjB49mpdffpkhQ4Zss5/33nuPww8/nO9///tccskl3HTTTVx++eXbHb979+6cc845AFx++eXccsstXHDBBdvFtWXLFp599lmmT5/OlVdeySOPPMItt9xCt27dmDlzJps2bWLkyJEcf/zxvPjii8ybN49XXnmF5cuXM2DAACZOnFjQaW8o16DMzKrZaaedOPPMM6smGqz09NNPc/rppwPJNBZPPPFEQfs75ZRTkMTgwYPZbbfdGDx4MG3atGHgwIEsWrQIgHvuuYcRI0YwfPhw5syZk7fJrkOHDlWTIh588MFV21Y3e/ZsPvKRjzB48GDuvPNO5szJ3znsU5/61Hb7euihh7j99tsZNmwYhx12GKtWreIf//gHjz32GOPHj6dt27bsueee9R6HsCFaRg1q80Z4ZSp7vvkOO+3QPm+RLR12Yk3vo5s5MDNrsAJqOk3poosuYsSIEZx99tk1lil0WovK6TLatGmzzdQZbdq0YcuWLSxcuJDrrruOmTNnsvPOOzNhwoS80260b9++6pg1TeEBSbPifffdx9ChQ5k8efI2o7Hniyt3XxHBT3/60+0GmZ0+fXqzT+PRMmpQO+8NXXZjU6eebN6hV95Hu/c/nJlzxdpNJQzWzMrBLrvswmc/+1luueWWqmVHHnkkd999N5DM5XTUUUcB20+jUV/vvvsuO+64I926dWP58uX1nmCw+vHXrl3LHnvswebNm+s93cYJJ5zAz3/+86qJGF9//XXee+89PvrRj3L33XezdetWli1bxowZM+q134ZoGQmqnlaue7/UIZhZGfjGN76xTW++n/zkJ9x2220MGTKEO+64g+uvvx5Ipoi/9tprGT58+HadJAoxdOhQhg8fzsCBA5k4cSIjR46s1/bVj3/11Vdz2GGHcdxxx9VrCg9IurUPGDCAESNGMGjQIL70pS+xZcsWPvnJT7LffvsxePBgzjvvvKoZe5tS+Uy38dbzsMu++dfvMwr2P5En56+ke+cOeYu037CCVf1P4ao/zmHusrVMOefwBsVdF0/NYdZwnm6j/BVzuo3yuAY1eFwyH1Q+/16Y/Nz/xOaLpxb1na23Nh6Vwsxas/JIUBVnQ8eu+ad8//Ok5o+nFvWdrbc2nsnXzFqzghKUpBOB64G2wM0RcU219XsBvwK6p2UmRcT0IsdaVPUZk8+jn5s1n4ho9t5iVhzFvmRUZ4KS1Ba4ATgOWALMlDQtInI76V8O3BMRP5c0AJgO9CtqpA1w5zy46/XKV72AZ6rWzV22fY+bnl060Ktrx22WLV61HljpBGXWDDp16sSqVavo0aOHk1SZiQhWrVpFp061D9hdH4XUoA4F5kfEAgBJdwOnArkJKoCd0ufdgLeKFmEjnHFA8oAPO0kAjL/pmYI7STRk9HMza5g+ffqwZMkSVqyo4ZqzZVqnTp3o06dP0fZXSILqDbyZ83oJcFi1MlcAD0m6ANgR+FhRojOzVqV9+/b079+/1GFYRhSSoPLVs6s3NI4HJkfE/0o6ArhD0qCI+GCbHUnnAucC7LXXXg2Jt8HabN1Ij4V/SF/1ynleu/YbugEUXN4jVpiZFUchCWoJ0DfndR+2b8L7AnAiQEQ8LakT0BPYZrjbiLgRuBGS+6AaGHODbOrSd5vXm3foVdB20bZ+5dtvcNOEmVkxFDKSxExgP0n9JXUATgOmVSvzT2A0gKSDgE6Av6nNzKzB6kxQEbEFOB94EJhL0ltvjqSrJI1Ni30DOEfSS8AUYEKUaoiKApy+f6kjMDOzuhR0H1R6T9P0asu+m/P8VaB+g0eVUGXPPjMzy65WOVismZllnxOUmZllUnmMxddK1TbwrAeSNbOWzgkqw2obeNYDyZpZS+cmPjMzyyQnKDMzyyQ38RXZtkMq1cxDIpmZ1c4JqgAL1sCkpwore3TvvozZu+5yHhLJzKx2TlB1OLp34WUXrEl+FpKgzMysdk5QdRizd+EJp9BalpmZ1c2dJMzMLJOcoMzMLJOcoMzMLJOcoMzMLJOcoMzMLJNaRi++fy+EP09i0IbNtGubP+eu2f1I3ukzupkDMzOzhir/BLXPqDqLdFq7GMAJysysjJR/gtr/xOQBzJ6/ku6dO2xXpN9zVzd3VE2utqk4moun/DCzplT+CapMVY7Zd/sbnTlz3/XbrCtknL7apuJoLp7yw8yakhNUiWzq0heAXy+E8YO2TTYep8/MzL34zMwso5ygzMwsk5ygzMwskwq6BiXpROB6oC1wc0Rck6fMZ4ErgABeiojTixhns+n1xlRW7DuuwdvXZ+6oStXLa2s3Ns+Zw8h9ezL6oN0aHIuZWTmrM0FJagvcABwHLAFmSpoWEa/mlNkP+DYwMiLekVS2fY93XfD7Bieo2uaOWr4e3t6Qf90rq6ov6UD7tesAnKDMrNUqpAZ1KDA/IhYASLobOBV4NafMOcANEfEOQES8XexAy0F95o6qdNIf4E+nbLus/YYVfG3OPsULzMysDBVyDao38GbO6yXpslz7A/tLelLSM2mT4HYknSvpOUnPrVjhrtRmZlazQhKU8iyLaq/bAfsBo4DxwM2Sum+3UcSNEVERERW9evWqb6xmZtaKFNLEtwTom/O6D/BWnjLPRMRmYKGkeSQJa2ZRoiyCTmsXFzzkUWOGRvKgtGZmxVFIgpoJ7CepP7AUOA2o3kPvPpKa02RJPUma/BYUM9DGWLP7kdsta79hBR02rsxbfsd35m637P1OPdm8Q+21Pg9Ka2ZWPHUmqIjYIul84EGSbua3RsQcSVcBz0XEtHTd8ZJeBbYCF0fEdn3TSuWdPqMLThoDHz6dOcfd1aDjNKTmdfr+DTqUmVmLV9B9UBExHZhebdl3c54H8PX0YfVwxgGljsDMLJs8WGwGtdm6kfYbkgpoj4V/2GZdISOdm5m1BE5QGbSpS1+ibfK8+nUvj3RuZq2Fx+IzM7NMcoKq5u19PlXqEMzMDCeo7TRmoFgzMyseJygzM8skd5KwBtu4ZSvTZi0tdRhmlnFtOnbZqSHbOUFZg/XdecdSh2Bm5aBNm7YN2qzYcZiZmRVDi6pBde7YjtXr36+z3PtbP2DXrp2aISIzM2uoFpWghvfdboaPvJ6cn3+Q2KzJN328p4M3s9aiRSWoLCh0Wo+6puWobfr4xavWAyudoMysRXOCKqJ803rkU8i0HDVNH99+wxq+NqdHg+IzMysn5ZOgOnWDdcvzr9u8EXbO823ezAqd1qMxEyKambUW5ZOg9juu5nWvTG2+OMzMrFm4m7mZmWVS+dSgDKh9rqhi8HxTZpYVTlBlpra5oorB802ZWVa0ygRV6A29lXxjr5lZ82uVCarQG3orlcuNvWZmLYk7SZiZWSY5QZmZWSYVlKAknShpnqT5kibVUm6cpJBUUbwQzcysNaozQUlqC9wAjAEGAOMlDchTritwIfD3YgdpZmatTyE1qEOB+RGxICLeB+4GTs1T7mrgh8DGIsZnZmatVCEJqjfwZs7rJemyKpKGA30j4o9FjM0yYOrzb9ZdyMysCRTSzVx5lkXVSqkN8GNgQp07ks4FzgXYa6+9CouwhaptWo66puJoTr97YSnjDu5b6jDMrBUqpAa1BMj9huoDvJXzuiswCHhU0iLgcGBavo4SEXFjRFREREWvXsUfBaFcrNn9SDZ2zT/6eqe1i+n2r6fyrjMza00KqUHNBPaT1B9YCpwGnF65MiLWAD0rX0t6FPhmRDxX3FBbjtqm5fBUHGZmiTprUBGxBTgfeBCYC9wTEXMkXSVpbFMHaGZmrVNBQx1FxHRgerVl362h7KjGh2V1WbAGJjVBS6C2dmPznDnbLLvqj3Pylh25b09PO29mTaZVjsVX7o7uXXeZ+lq+Ht7eANABVq/dZt3cZWu3K9++bdJ3xgnKzJqKE1QZGrN38mgK7TesYFX/U6pej7/pGaacc/h25WqqVZmZFYsTVAHqMz2Hp+YwMysOJ6gC1Gd6jnKfmqPN1o3VZuptvbcDmFlptYwE1akbrFtee5nNG2HnJmoXa0E2dfFNuWaWDS0jQe13XN1lXpna9HGYmVnReD4oq9V/9n+v1CGYWSvlBGW1OnPf9aUOwcxaKScoMzPLJCcoMzPLJCcoMzPLpJbRi6+FKZe5oszMmpITVMas2f3IGtd1WrsYwAnKzFoFJ6iM8VxRZmYJX4MyM7NMaj01qEKGQyqCzu8tg877NPlxzMxautaToAoZDqkIOrx1a8Ejn9fFI6ObWWvWehJUMxm4507QpWdR9lXuI6ObmTWGr0GZmVkmOUGZmVkmOUGZmVkmOUGZmVkmOUGZmVkmFZSgJJ0oaZ6k+ZIm5Vn/dUmvSnpZ0l8keW51MzMDoE2nLjs3aLu6CkhqC9wAjAEGAOMlDahW7EWgIiKGAFOBHzYkGDMza3nadOzcNAkKOBSYHxELIuJ94G7g1NwCETEjIiqnXn0G6NOQYMzMzCoVcqNub+DNnNdLgMNqKf8F4IF8KySdC5wLsNdeexUYouWqbSqOYvB0HmaWFYUkKOVZFnkLSv8JVABH51sfETcCNwJUVFTk3UfZK3TMv80bYef6XaqrbSqOYvB0HmaWJYUkqCVA35zXfYC3qheS9DHgMuDoiNhUnPDKUKFj/r0ytd67rm0qjmKob81s8ar1XPXHOU0UjZm1doUkqJnAfpL6A0uB04DTcwtIGg78EjgxIt4uepSWOSP37Ql4rEAz+9CKtZtYua44g2VDAQkqIrZIOh94EGgL3BoRcyRdBTwXEdOAa4EuwG8lAfwzIsYWLUrLnNEH7cbog3YrdRhmVgZG/aph2xU0mnlETAemV1v23ZznH2vY4a02nTu2q3HqDk/FYWYtnafbyLDhfbvXuM5TcZhZS+ehjszMLJOcoMzMLJOcoMzMrEl9sGn9Ow3ZzgnKzMya1Acb1zlBmZlZy+FefKVS25BIDRgGycyspXGCKpXahkRqwDBIZmYtjZv4zMwsk1yDsm1Un85jXY8hrOp/SgkjMrPWygnKqlSfzqPT2sVo6+YSRWNmrZ0TVJmqbZy+hlq9y0dgl49UvR70yv/Q5mhb04EAABAvSURBVIPNrFi7MW/5jVu20nfnHYsag5lZJSeoMlXbOH1FM789bIWxw3rnXT1t1tKmj8HMWi13kjAzs0xyDSqLfI+UmZkTVCb5HikzMzfxmZlZNjlBmZlZJjlBmZlZJjlBmZlZJjlBmZlZJrkXX7mprQt6Mbgbu5llREEJStKJwPVAW+DmiLim2vqOwO3AwcAq4HMRsai4oRpQexf0YnA3djPLiDqb+CS1BW4AxgADgPGSBlQr9gXgnYj4D+DHwA+KHaiZmbUuhdSgDgXmR8QCAEl3A6cCr+aUORW4In0+Ffh/khQRUcRYLWO67tC+xoFkzcyqfPDB1oZsVkiC6g28mfN6CXBYTWUiYoukNUAPYGVuIUnnAucC7LXXXg2J15pa7jWurrtD2441Fj3mgF2bKSgzK2cfbFr3bkO2KyRBKc+y6jWjQsoQETcCNwJUVFS4dpVFude4Bo8rXRxm1uoV0s18CdA353Uf4K2aykhqB3QD/l2MAM3MrHUqJEHNBPaT1F9SB+A0YFq1MtOAs9Ln44C/+vqTmZk1Rp1NfOk1pfOBB0m6md8aEXMkXQU8FxHTgFuAOyTNJ6k5ndaUQZuZWctX0H1QETEdmF5t2Xdznm8EPlPc0MzMrDXzUEdmZpZJTlBmZpZJKlVfBklrgXklOXjj9aTaPV5lppzjd+yl4dhLp5zjr4x974joVd+NSzlY7LyIqCjh8RtM0nPlGjuUd/yOvTQce+mUc/yNjd1NfGZmlklOUGZmlkmlTFA3lvDYjVXOsUN5x+/YS8Oxl045x9+o2EvWScLMzKw2buIzM7NMcoIyM7NMKkmCknSipHmS5kuaVIoYCiWpr6QZkuZKmiPpa+nyXSQ9LOkf6c+dSx1rTSS1lfSipD+mr/tL+nsa+2/SQYAzR1J3SVMlvZae/yPK5bxL+q/08zJb0hRJnbJ83iXdKultSbNzluU910r8JP37fVnSiNJFXmPs16afm5cl3Supe866b6exz5N0Qmmiroplu9hz1n1TUkjqmb7O/HlPl1+Qnts5kn6Ys7z+5z0imvVBMuDsG8A+QAfgJWBAc8dRj3j3AEakz7sCrwMDgB8Ck9Llk4AflDrWWt7D14G7gD+mr+8BTkuf/wI4r9Qx1hD3r4Avps87AN3L4byTTOC5ENgh53xPyPJ5Bz4KjABm5yzLe66BjwMPkMwDdzjw9wzGfjzQLn3+g5zYB6TfOR2B/ul3UdssxZ4u70syQPdioGcZnfdjgEeAjunrXRtz3ktRg6qaQj4i3gcqp5DPpIhYFhEvpM/XAnNJvoBOJfkCJf35idJEWDtJfYCTgJvT1wKOBaamRTIZu6SdSP4AbgGIiPcjYjVlct5JboLfIZ0frTOwjAyf94h4jO3ncKvpXJ8K3B6JZ4DukvZonki3ly/2iHgoIrakL58hmccOktjvjohNEbEQmE/ynVQSNZx3gB8Dl7DtxK+ZP+/AecA1EbEpLfN2urxB570UCSrfFPK9SxBHvUnqBwwH/g7sFhHLIEliQFbnP/8/kg/6B+nrHsDqnD/erJ7/fYAVwG1p8+TNknakDM57RCwFrgP+SZKY1gDPUx7nPVdN57rc/oYnktQ8oAxilzQWWBoRL1VblfnYgf2Bj6RN2X+TdEi6vEGxlyJBFTQ9fNZI6gL8DrgoIt4tdTyFkHQy8HZEPJ+7OE/RLJ7/diTNBz+PiOHAeyTNTJmXXqs5laQpY09gR2BMnqJZPO+FKJfPEJIuA7YAd1YuylMsM7FL6gxcBnw33+o8yzITe6odsDNJE+TFwD1pq02DYi9FgipkCvlMkdSeJDndGRG/Txcvr6xepz/frmn7EhoJjJW0iKQp9ViSGlX3tOkJsnv+lwBLIuLv6eupJAmrHM77x4CFEbEiIjYDvweOpDzOe66aznVZ/A1LOgs4GTgj0gshZD/2fUn+sXkp/bvtA7wgaXeyHzskMf4+bYZ8lqTlpicNjL0UCaqQKeQzI83+twBzI+JHOatyp7k/C7i/uWOrS0R8OyL6REQ/kvP814g4A5gBjEuLZTX2fwFvSjogXTQaeJUyOO8kTXuHS+qcfn4qY8/8ea+mpnM9DTgz7VV2OLCmsikwKySdCHwLGBsR63NWTQNOk9RRUn9gP+DZUsSYT0S8EhG7RkS/9O92CUknrX9RBucduI/kH2Ek7U/SuWklDT3vJer98XGS3nBvAJeVIoZ6xHoUSVX0ZWBW+vg4ybWcvwD/SH/uUupY63gfo/iwF98+6YdjPvBb0h43WXsAw4Dn0nN/H0nTQVmcd+BK4DVgNnAHSe+lzJ53YArJ9bLNJF+KX6jpXJM019yQ/v2+AlRkMPb5JNc8Kv9mf5FT/rI09nnAmKzFXm39Ij7sxVcO570D8Ov0c/8CcGxjzruHOjIzs0zySBJmZpZJTlBmZpZJTlBmZpZJTlBmZpZJTlBmZpZJTlBmZpZJTlDWYiiZnuMrOa/3lDS1tm3qse8rJH0zfX6VpI8Vab97KJ0GpSlIWlePso8oo9OXWOvkBGUtSXegKkFFxFsRMa6W8g0SEd+NiEeKtLuvAzcVaV+NdQc558+s1JygrCW5BthX0qx0wrp+lZOpSZog6T5Jf5C0UNL5kr6ejpT+jKRd0nL7SvqzpOclPS7pwOoHkTRZ0rj0+SJJV0p6QdIrleUl7ZhO6DYzPUZNU8p8Gvhzus10SUPS5y9K+m76/GpJX0yfX5zu82VJV+bE9J+Snk3f+y8lta0Wc09JT0s6Ka21PZaWnS3pI2mxacD4Bp57s6JzgrKWZBLwRkQMi4iL86wfBJxOMg/N94H1kYyU/jRwZlrmRuCCiDgY+CbwswKOuzIiRgA/T7eBZFiXv0bEISSTuF2bThdSJR2T7J1I584BHiOZqmAnkhG4R6bLjwIel3Q8yRhmh5IMA3WwpI9KOgj4HDAyIoYBW4Ezco6zG/An4LsR8af0HDyYlh1KMhQQEfEO0FFSjwLes1mTa1d3EbMWY0Ykk06ulbQG+EO6/BVgSDqlypHAb5MxXoFkDL26VI5w/zzwqfT58SQjyVcmrE7AXiQTXlbag2TOq0qPAxeSzMb7J+C4dPqFfhExT9I56X5fTMt3IUlYQ4CDgZlp3Dvw4cjj7UnG0ftqRPwtXTYTuDUdpf++iJiVE8PbJFOErCrgfZs1KScoa0025Tz/IOf1ByR/C21IJhUc1sD9buXDvykBn46IebVst4EkcVWaCVQAC4CHSaYpOIck8VXu838i4pe5O5F0AfCriPh2nmNsSbc/AfgbJDOhSvooyUzLd0i6NiJuT8t3SuMyKzk38VlLshbo2tCNI5mIcqGkz0Ay1YqkoQ3c3YPABel0G0ganqfM60C/nOO/TzIC92dJpil/nKTJ8PGcfU5Ma3pI6i1pV5Ia0rj0OZJ2kbR35W5JZpQ9UNKkdP3eJBNZ3kQylcyIyvcL7E4ygrZZyTlBWYsREauAJ9ML/9c2cDdnAF+Q9BIwh2Rm3Ia4mqR57eW0o8bVeeJ9D3hD0n/kLH4cWB7JHEaPk0zs9nha/iHgLuBpSa+QTOLYNSJeBS4HHpL0Mknta4+c42wlmQ/smLQb/ihglqQXSTppXJ8WPRh4Jj6clt6spDzdhlkJSfokcHBEXJ6BWK4HpkXEX0odixkUoQYl6ZOSIrc7btq99/R67OOpOtYXfLOhWTmJiHvJTpPabCcny5JG16Ak3UPSnPCXiLgiXTYK+GZEnFzHtm3T5oe6jrEuIro0KlAzMysrjapBpRdrR5JM9XtazqprSO7nmCXpv6ptM0rSDEl3kXTvraoh1XIDYeW2VTcbNiZuMzPLvsZ2M/8E8OeIeF3SvyWNiIgXSG6YrK0GdSgwKCIWVlteeQPh99M74TtXrkhvNpwGXB4RDzcybjMzy7jGXoMaD9ydPr+bwodJeTZPcoLkPpCzJV0BDE5vqoQPbza8xMnJzKx1aHCCSodDORa4WdIi4GLgc5X3fdThvXwLI+Ix4KPAUpIbCCuHn8m92dDMzFqBxtSgxgG3R8TeEdEvIvqSDNFyFA28YbKmGwjJc7OhmZm1bI1JUOOBe6st+x3JdaSXgS2SXqreSaIOo8h/A2G+mw3NzKwF8426ZmaWSR7qyMzMMskJyszMMskJyszMMskJyszMMskJyszMMskJyszMMqloCUrSrZLeTidnq1z2GUlzJH0gqaJYx2pOkg5IB6+tfLwr6aJSx1UISZ0kPZvejzZH0pWljqlQ+T5P5aSc43fspVHOsUPTxF/MGtRk4MRqy2YDnwIeK+JxmlVEzIuIYRExjGTG0fVsf4NyVm0Cjo2IocAw4ERJh5c4pkJNZvvPUzmZTPnGPxnHXgqTKd/YoQniL1qCSsfR+3e1ZXMjYl6xjpEBo4E3ImJxqQMpRCQqJ3tsnz7K4s7sfJ+nclLO8Tv20ijn2KFp4vc1qPo5DZhS6iDqQ1JbSbOAt4GHI+LvpY7JzKwQTlAFktQBGAv8ttSx1EdEbE2bJ/sAh0oaVOqYzMwK4QRVuDHACxGxvNSBNERErAYepbzbuM2sFXGCKtx4yq95r5ek7unzHYCPAa+VNiozs8IUs5v5FOBp4ABJSyR9QdInJS0BjgD+JOnBYh2vOUnqDBwH/L7UsdTTHsAMSS+TzFb8cET8scQxFSTf56nUMdVHOcfv2EujnGOHponf022YmVkmuYnPzMwyyQnKzMwyyQnKzMwyqanH4rtW0muSXpZ0b2WPsnIjqbukqel7mSvpiFLHVAhJfSXNSGOeI+lrpY6pPiSdKGmepPmSJpU6nvoo59ihvON37KXRJLFHRFEewEeBEcDsnGXHA+3S5z8AflCs4zXnA/gV8MX0eQege6ljKjDuPYAR6fOuwOvAgFLHVWDsbYE3gH3Sc/6SY3f8jj2bj6aKvanH4nsoIrakL58hGc2grEjaiST53gIQEe9HctNr5kXEsoh4IX2+FpgL9C5tVAU7FJgfEQsi4n3gbuDUEsdUqHKOHco7fsdeGk0Se3Neg5oIPNCMxyuWfYAVwG2SXpR0s6QdSx1UfUnqBwwHymUsvt7Amzmvl1A+ybWcY4fyjt+xl0aTxN4sCUrSZcAW4M7mOF6RtSNpuvx5RAwH3gPKrW24C/A74KKIeLfU8RRIeZaVy0175Rw7lHf8jr00miT2Jk9Qks4CTgbOiLSxsswsAZbEh6OATyVJWGVBUnuS5HRnRJTTSBhLgL45r/sAb5Uolvoq59ihvON37KXRJLE3aYKSdCLwLWBsRKxvymM1lYj4F/CmpAPSRaOBV0sYUsEkieTa2dyI+FGp46mnmcB+kvqnI8mfBkwrcUyFKufYobzjd+yl0SSxt2t0WKl0HKZRQM90/L3vAd8GOgIPJ9+VPBMRXy7WMZvRBcCd6YlfAJxd4ngKNRL4PPBKOicUwKURMb2EMRUkIrZIOh94kKSH0K0RMafEYRWknGOH8o7fsZdGU8XusfjMzCyTPJKEmZllkhOUmZllkhOUmZllkhOUmZllkhOUmZllkhOUmZllkhOUmZll0v8PF9DvOdWKLbcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEOCAYAAADc94MzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwU1b3//9ebXQVBWdwwgrluIKuouERR3BeMiUlcEkWMJt6oWVFcrnG55muiN7lJfmZxRY1KvCQiUYxLXDCuoKKCiqKgjhocCSDgxsDn90fVjE3TM9NAN1Uw7+fj0Y+p5fSpz9R0z6fPqepzFBGYmZnlTausAzAzMyvFCcrMzHLJCcrMzHLJCcrMzHLJCcrMzHLJCcrMzHLJCcqqRtIcSQfkII4TJN1XgXpGSvpnJWJqySQNk1TTxP6Q9B9rM6Y8kDRW0n9nHUeeOEFVQPqP+GNJiyXNlXSDpI6rWVeTb94y6/ihpH9JWijpeknt16S+apDUU9JfJH2QxvmipJHVOFZE3BIRB1Wj7pZoXflHKmkLSRMlvZsmvV4lyhwg6VlJSyS9Lenraz9Sa4wTVOUcGREdgcHArsAFq1qBpDZrGoSkg4ExwHCgF7AtcPGa1lsFNwNvA9sAXYETgbmrU1Elztv6xOejwXLg78BXS+2U1Ae4FTgf6AwMBJ5Za9FZs5ygKiwi3gHuAXYGkHSypJclLZL0hqTv1Jetby1JOkfSv4Db0udumbbGFkvaUtJHkroWPG8XSbWS2pYI4STguoiYERHzgUuBkeXGL+n/ClpfkyX1Ldg3VtLvJN2TxvaYpM0l/a+k+ZJekTSozEPtCoyNiCURURcRz0XEPYXnpSiuhu5CSRdJGi/pT5I+BM5LW7CbFpQflLbO2hZ2zUn6g6Qri+q+U9KP0uUxkl5P/14vSTq63HPXnKbqltRa0v+kMc+WdEb6qb9Nur93+vdYJOkBSVdJ+lO6r1da9hRJbwEPptuHSnpc0gJJz0saVnC8RutL95d8HUg6DTgBODt9Dfwt3b5l2iKuTeM/q6CuDdLXznxJL5H87ZtzWPp++UDSFZJaSWov6d+S+hXU3SP923cvriAi5kbE74ApjRzjAuCPEXFP+hqcFxGvlyoo6RFJX02X907P92Hp+gGSphWUHZW+5+dLulfSNgX7dpR0f/p7zFQjLTZJnSQ9JOk3ShyWvmYWSXpH0k/KOIfrPCeoCpO0NXAY8Fy66X3gCGBj4GTgV5IGFzxlc2BTkpbEicChwLsR0TF9vAs8DBS+kL8JjIuIpSVC6As8X7D+PLCZ0gQn6S5JY5r4Fe4BtgN6AM8CtxTt/zrJG7sb8CnwRFquGzAe+GUTdRd6ErhK0rGSvlDmcwodlR6vC3BFGkfhJ+XjgfElztGtwDckCUDSJsBBwLh0/+vAl0g+UV8M/EnSFqsRXylN1X0qyd9+IEkr/Msl4n6apLV5EfCtEvXvC+wEHCxpK+Bu4L9JXl8/Af5S8I+8ufpKvg4i4up0+Rfp6/NISa2Av5G81rYiab3/QElrHuCnwBfTx8EkH6KaczQwJD0XRwGjIuJTkr/TNwvKHQc8EBG1ZdRZbCiAku7l99IPPJs2UvYRYFi6vA/wBsn5rl9/JK3ry8B5wFeA7sCjJB88kbQRcD/Jue+Rxv47FXwITMt1Bf4BPBYRZ0UyHt11wHciohPJh98HV+P3XfdEhB9r+ADmAIuBBcCbwO+ADRopOwH4fro8DPgM6FCwfxhQU/Scb5C8WAFaA/8Cdmuk/teBQwrW2wIB9FqN36tL+tzO6fpY4JqC/WcCLxes9wMWFJ2XAxqpexPgcmAGsAyYBuzaxDloqIvkH+rkov3fBh5Ml0XSfbhPuj4S+GfBvrcK9p1a/7xG4pwGHFVcT4VeN4V1P0jyD6h+3wHpuW8DfAGoAzYs2P8n4E/pcq+07LYF+88Bbi463r0kyaHJ+sp8Hfx3wf7dgbeKnnMucEO6/EbRa/K04r9v0XOjqPx/Av8oONbbQKt0fSrw9WbOcxtKvAdI3ntzgO2BjsBfgFsaqWM48EK6/Pf09fZkuv4I8JV0+R7glILntQI+IvkA+g3g0aJ6/wj8tOC8Xg9MB0YXlXsL+A6wcaVef+vCwy2oyvlyRHSJiG0i4j8j4mMASYdKejJt0i8gaV11K3hebUR80kzddwJ9JG0LHAgsjIinGym7mKS1Vq9+eVFzv0DazXR52g31Icmbl6J4C68TfVxivaybQyJifkSMiYi+wGYk/6wn1LdsyvB20fp4YA9JW5J8og2ST6/Fxw2ST+HHpZuOp6CVKOlESdPSbrEFJJ9WuxXXU0xJ12F9t+x5jZRpqu4ti36nwuUtgX9HxEeN7C+1bRvga/XHSo+3N7BFc/WV+TootA1Jt3Thsc4j+buW+t3ebKSexn6XN9M6iIingCXAvpJ2BP4DmFhGfaV8TJJEX42IxcDPSN6fpTwBbC9pM5JW7k3A1pK6AbsBk9Ny2wC/LjgP/yb5ULRVum/3ovN0AkkvSr3DgQ2APxQd/6tpbG+m3Y17rObvvE7xxdQqUnL33F9Iuu7ujIilkiaQvGDrFQ8nv9Lw8hHxiaTbSV7MO5LcYNCYGcAA4PZ0fQAwNyLmlRHy8STdKQeQ/FPqDMwvirfiIuIDJdeFTiLpjloCbFi/X1Jrku6SFZ5WVMcCJbeSf52km+u2NBmVchtwn6TLST6RH50eZxvgGpJPy09ExLL02kKzv39EfBf4bmP7y6j7PaBnwVO2Llh+D9hU0oYFSaVwf0MYBctvk7SgTm0klqbqa+51UHxe3wZmR8R2JWKqj39rktcmJC245hSXf7dg340k3Xz/IunGbe4DXmNeoMT7rZSI+EjSM8D3gekR8Zmkx4EfAa9HxAdp0beByyKiuGu8/rw/EhEHNnGoa0h6FyZJOiQilqTHnwIcpeS68xkk7+9Sr4H1iltQ1dUOaA/UAnWSDiW53tGUuUBXSZ2Ltt9E0sU0gqQ7pjE3AadI6pNeX7mApOugHJ1IrivNI0kQPyvzeatM0s8l7SypjaROwOnArDSRvgp0kHR4+oa8gOQ8NudWkg8DX02XS4qI50j+JtcC90bEgnTXRiT/sGrTGE8mvdmlApqr+3bg+5K2ktSFpIuuPt43SbqyLpLULv30fGQzx/sTcKSkg9MWUQclN5/0LKO+5l4Hc0nuDq33NPChkpt9NkiPt7Ok+pshbgfOlbSJpJ4kXcPNGZ2W35okKfy5YN/NJB8qvknyem+UpA58/tppn67XuwE4WdK2kjYkOed3NVHdIyTJ4ZF0/eGidUhaPufq85tKOkv6WrrvLpJW2LeU3LzTVtKuknYqOs4ZwEzgrvR8tlPyXb7OkVxT/ZCkW3y95wRVRRGxCDiL5A06n+STaZPdERHxCskn/DfSboD6ro3HSG6bfTYi5jTx/L8DvwAeIukaeZPkIjUASu7AK9kFRfJmfxN4B3iJ5EaGatkQuIPkut0bJN0fIwAiYiHJdYdr01iWAOV8N2wiyYX9uRHxfDNlbyNpITQksoh4Cfgfku6cuSTX1B4r+zdqQhl1XwPcR/Kp/jlgEsl1ovp/RCcAe5Akjf8m+Yf9aRPHe5ukFXQeSVJ8GxjN5+/5pupr7nVwHUmX8wJJEyJiGUmCGwjMBj4g+dvVf8i6OK1vdvo7NtUDUO9Oklu+p5Hc7HFdwe9WQ3LjRslu3CIfk3R7A7ySrtfXc336uz6Vxvcpyfu1MY+QJO/JjawTEXcAPwfGpd2j00lufqn/f3AQcCxJi/BfadkVPnylLf/TSP5mdwIdSG5imZPW+V1WvFFkvaXGe0EsbyQ9CNwaEddmHYtVV9ra/kNEbNPI/j8Dr0TET0vtX43jVbS+apN0Pcndrqv8fUNbd7gFtY5Iu0sGs2JXh60n0q6cw9Iuz61IWr13FOzfVdIXlXwf6BCS1tGENTheRetbm5SMCPEVClpVtn5ygloHSLoReAD4QdpNYOsfkXSFzSfp4nsZuLBg/+Yk1zwWA78BTk+vpa2uSte3Vki6lKTb7IqImJ11PFZd7uIzM7NccgvKzMxyKbPvQXXr1i169eqV1eHNzCwnnnnmmQ8iYqXxFDNLUL169WLq1KlZHd7MzHJCUsnRRdzFZ2ZmueQEZWZmueQEZWZmudTsNaj0G9tHAO9HxErjkkkS8GuSkXY/AkZGxLOVDtTMWralS5dSU1PDJ5+s7tiwlrUOHTrQs2dP2rYtNdfqysq5SWIs8P/R+KCMh5KMf7YdycjQv09/mplVTE1NDZ06daJXr16o7FlZLC8ignnz5lFTU0Pv3r3Lek6zCSoiJqdDizTmKOCmdIDDJyV1kbRFRLzXVL0fv/cKM362d1lBVtvi7Y5m96/9OOswzKwJn3zyiZPTOkwSXbt2pba2/MmPK3ENaitWnFysJt22EkmnSZoqaWpeRrDY+rPX6fjaHc0XNLPMOTmt21b171eJ70GVOmLJ7BMRVwNXAwwZMiT6nvfPChx+zeSlFWdmZiuqRAuqhhVnduzJirNfmpmtFyTxrW99q2G9rq6O7t27c8QRRzT5vKlTp3LWWU1NNQULFizgd7/7XVlx7LnnnmWVa86cOXPYeedKzclZeZVIUBOBE5UYCixs7vqTmdm6aKONNmL69Ol8/HEy7+H999/PVluVvKKxgiFDhvCb3/ymyTKrkqAef/zxssqt65pNUJJuI5kFdAdJNZJOkfRdSd9Ni0wimRF1FsmsoP9ZtWjNzDJ26KGHcvfddwNw2223cdxxxzXse/rpp9lzzz0ZNGgQe+65JzNnzgTg4YcfbmhlXXTRRYwaNYphw4ax7bbbNiSuMWPG8PrrrzNw4EBGjx7N4sWLGT58OIMHD6Zfv37ceeedDcfp2LFjQ73Dhg3jmGOOYccdd+SEE06g/vr+M888w7777ssuu+zCwQcfzHvvvdewfcCAAeyxxx5cddVVVT5ba6acu/iOa2Z/AN+rWERmZs24+G8zeOndDytaZ58tN+anR/Ztttyxxx7LJZdcwhFHHMELL7zAqFGjePTRZOb5HXfckcmTJ9OmTRseeOABzjvvPP7yl7+sVMcrr7zCQw89xKJFi9hhhx04/fTTufzyy5k+fTrTpk0Dku7DO+64g4033pgPPviAoUOHMmLEiJVuNHjuueeYMWMGW265JXvttRePPfYYu+++O2eeeSZ33nkn3bt3589//jPnn38+119/PSeffDK//e1v2XfffRk9enQFzlz1ZDZYrJnZuqh///7MmTOH2267jcMOO2yFfQsXLuSkk07itddeQxJLly4tWcfhhx9O+/btad++PT169GDu3LkrlYkIzjvvPCZPnkyrVq145513mDt3LptvvvkK5XbbbTd69uwJwMCBA5kzZw5dunRh+vTpHHjggQAsW7aMLbbYgoULF7JgwQL23XdfAL71rW9xzz33rPE5qRYnKDNb55TT0qmmESNG8JOf/ISHH36YefPmNWz/r//6L/bbbz/uuOMO5syZw7Bhw0o+v3379g3LrVu3pq6ubqUyt9xyC7W1tTzzzDO0bduWXr16lRxFo1RdEUHfvn154oknVii7YMGCdepWfScooG5ZMHHaO1mH0aDTBm3Zb4ceWYdhZo0YNWoUnTt3pl+/fjz88MMN2xcuXNhw08TYsWNXqc5OnTqxaNGiFerq0aMHbdu25aGHHuLNN0vOSFHSDjvsQG1tLU888QR77LEHS5cu5dVXX6Vv37507tyZf/7zn+y9997ccsstqxTj2uYElereqUPWITSoXeSxxszyrGfPnnz/+99fafvZZ5/NSSedxC9/+Uv233//Vaqza9eu7LXXXuy8884ceuihnHPOORx55JEMGTKEgQMHsuOOO5ZdV7t27Rg/fjxnnXUWCxcupK6ujh/84Af07duXG264gVGjRrHhhhty8MEHr1KMa5uyGtFhyJAhkYcJC2f8bG/qlgUfHT8x61Aa1C76hBEDm7911awlefnll9lpp52yDsPWUKm/o6RnImJIcVlPt2FmZrnkBGVmZrnkBGVmZrnkBGVmZrnkBGVmZrnkBGVmZrnkBGVmViZJ/PjHn8++feWVV3LRRRc1+ZwJEybw0ksvVSWeb3/7283WXe7x//CHP3DTTTdVJK6RI0cyfvz4Na7HCcrMrEzt27fnr3/9Kx988EHZz6lmgrr22mvp06dPRY7/3e9+lxNPPLFSoVWEE5SZWZnatGnDaaedxq9+9auV9r355psMHz6c/v37M3z4cN566y0ef/xxJk6cyOjRoxk4cCCvv/76Cs8ZOXIkp59+Ovvttx/bbrstjzzyCKNGjWKnnXZi5MiRDeVOP/10hgwZQt++ffnpT3/asH3YsGHUD3jQsWNHzj//fAYMGMDQoUOZO3duyeNfc8017LrrrgwYMICvfvWrfPTRR0AyDciVV17ZUO8555zDbrvtxvbbb98wWvuyZcsYPXo0u+66K/379+ePf/wjkAxse8YZZ9CnTx8OP/xw3n///cqc74rUYma2Nt0zBv71YmXr3LwfHHp5s8W+973v0b9/f84+++wVtp9xxhmceOKJnHTSSVx//fWcddZZTJgwgREjRnDEEUdwzDHHlKxv/vz5PPjgg0ycOJEjjzySxx57jGuvvZZdd92VadOmMXDgQC677DI23XRTli1bxvDhw3nhhRfo37//CvUsWbKEoUOHctlll3H22WdzzTXXcMEFF6x0/C5dunDqqacCcMEFF3Dddddx5plnrhRXXV0dTz/9NJMmTeLiiy/mgQce4LrrrqNz585MmTKFTz/9lL322ouDDjqI5557jpkzZ/Liiy8yd+5c+vTpw6hRo8o67U1xC8rMbBVsvPHGnHjiiSvNkPvEE09w/PHHA8k0Fv/85z/Lqu/II49EEv369WOzzTajX79+tGrVir59+zJnzhwAbr/9dgYPHsygQYOYMWNGyS67du3aNUyKuMsuuzQ8t9j06dP50pe+RL9+/bjllluYMWNGyXJf+cpXVqrrvvvu46abbmLgwIHsvvvuzJs3j9dee43Jkydz3HHH0bp1a7bccstVHoewMW5Bmdm6p4yWTjX94Ac/YPDgwZx88smNlil3Wov66TJatWq1wtQZrVq1oq6ujtmzZ3PllVcyZcoUNtlkE0aOHFly2o22bds2HLOxKTwg6VacMGECAwYMYOzYsSuMxl4qrsK6IoLf/va3Kw0yO2nSpKpM4+EWFLB02fKsQzCzdcimm27K17/+da677rqGbXvuuSfjxo0Dkrmc9t57b2DlaTRW1YcffshGG21E586dmTt37ipPMFh8/EWLFrHFFluwdOnSVZ5u4+CDD+b3v/99w0SMr776KkuWLGGfffZh3LhxLFu2jPfee4+HHnpoleptjBMUULc8mxHdzWzd9eMf/3iFu/l+85vfcMMNN9C/f39uvvlmfv3rXwPJFPFXXHEFgwYNWukmiXIMGDCAQYMG0bdvX0aNGsVee+21Ss8vPv6ll17K7rvvzoEHHrhKU3hAclt7nz59GDx4MDvvvDPf+c53qKur4+ijj2a77bajX79+nH766Q0z9q4pT7fxs71Z9EkdcdLdWYfSwNNtmK3M022sHzzdhpmZrfOcoMzMLJfKuotP0iHAr4HWwLURcXnR/i8ANwJd0jJjImJShWOtqkvuKn2rZRb6b9XZXXxmJUREVe4Ws7VjVS8pNZugJLUGrgIOBGqAKZImRkThjfgXALdHxO8l9QEmAb1WKZK14Ff3v8qv//HaCtvGtUtun3z5vZXvsunWsR3dO7VfaXs1vTnvI5bW+a5Cs2IdOnRg3rx5dO3a1UlqHRQRzJs3jw4dOpT9nHJaULsBsyLiDQBJ44CjgMIEFcDG6XJn4N2yI1iLfnjg9vzwwO1X3HjD73hy9jxuO3VoNkEVueSuGU5QZiX07NmTmpoaamtrsw7FVlOHDh3o2bNn2eXLSVBbAW8XrNcAuxeVuQi4T9KZwEbAAWVHYGZWhrZt29K7d++sw7C1qJwEVaotXdyReBwwNiL+R9IewM2Sdo6IFZoCkk4DTgP4whe+sDrxVk3X2X/LOgQA2n7cGcL3rpiZlZOgaoCtC9Z7snIX3inAIQAR8YSkDkA3YIUhbSPiauBqSL4HtZoxV8XSDbpnHQIA0Rq09NOswzAzy1w5H9WnANtJ6i2pHXAsMLGozFvAcABJOwEdAHcUm5nZams2QUVEHXAGcC/wMsndejMkXSJpRFrsx8Cpkp4HbgNGRlZDVKyG7u2XZR2CmZkVKet7UOl3miYVbbuwYPklYNUGiMqRHu3reCvrIMzMbAW+Gm9mZrnk+aByKAgmTnsn6zAA6LRBW/bboUfWYZhZC+QElUNtW7Wie6fyv21dTbWLVp4YzcxsbXAXn5mZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5ZITlJmZ5VJZCUrSIZJmSpolaUwjZb4u6SVJMyTdWtkwq6vdZ/OzDsHMzIq0aa6ApNbAVcCBQA0wRdLEiHipoMx2wLnAXhExX1KPagVcDe2doMzMcqecFtRuwKyIeCMiPgPGAUcVlTkVuCoi5gNExPuVDdPMzFqachLUVsDbBes16bZC2wPbS3pM0pOSDilVkaTTJE2VNLW2tnb1IjYzsxahnASlEtuiaL0NsB0wDDgOuFZSl5WeFHF1RAyJiCHdu3df1VjNzKwFafYaFEmLaeuC9Z7AuyXKPBkRS4HZkmaSJKwpFYlyLeg19dKsQwDgZx/D5Fa7A/tkHYqZWabKSVBTgO0k9QbeAY4Fji8qM4Gk5TRWUjeSLr83KhloRTz0/+CRy0vu2mj+yytt+6xDN5ZusHZber2XvwmxnMVr9ahmZvnTbIKKiDpJZwD3Aq2B6yNihqRLgKkRMTHdd5Ckl4BlwOiImFfNwFfLfucmj2IXdWbGgfm4M375Py6FWJ51GGZmmSunBUVETAImFW27sGA5gB+lDzMzszVWVoKyta/r7L9lHUKirj1J762Z2drlBJVHarXWr3016oO3mTjtnayjMLP1WKv2HTcutd0JKoc+Ww4dsg4i1XWjDtApL9GY2XqpVavWJTev7Tjy6N0tDsw6hBUs9T0SZmZOUADvbnlw1iGYmVkRJygzM8slJygzM8sl3ySRU2MezzqCxP7dOzC0d9ZRmFlL5ASVoVtmwq2vrrhtXLvk54slxuHosQFstmH146r3xkLQsvYMXXuHNDNr4ASVoRN2SB6Fek1NktPdR2YTU6Exj5MMXGVmlgFfgzIzs1xygjIzs1xygjIzs1xygjIzs1xygsqhHhtkHYGZWfacoHJobd5KbmaWV05QZmaWS05QZmaWS/6iLrBBu9Ys+OizrMMAoG7ZcpZH1lGYmWXPCQrou+XG0LFb1mEkZrVl4cdLs47CzCxz7uIzM7NccoIyM7NcKitBSTpE0kxJsySNaaLcMZJC0pDKhWhmZi1RswlKUmvgKuBQoA9wnKQ+Jcp1As4Cnqp0kGZm1vKU04LaDZgVEW9ExGfAOOCoEuUuBX4BfFLB+MzMrIUqJ0FtBbxdsF6TbmsgaRCwdUTcVcHYzMysBSsnQanEtoZv6khqBfwK+HGzFUmnSZoqaWptbW35UZqZWYtTToKqAbYuWO8JvFuw3gnYGXhY0hxgKDCx1I0SEXF1RAyJiCHdu3df/ajNzGy9V06CmgJsJ6m3pHbAscDE+p0RsTAiukVEr4joBTwJjIiIqVWJ2MzMWoRmE1RE1AFnAPcCLwO3R8QMSZdIGlHtAM3MrGUqa6ijiJgETCradmEjZYeteVhmZtbSeSQJMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLpbLG4rO1a6Mlb9Fr6qVZh8HPPoa/LR8K7JN1KGa2HmvVoeMmpbY7QeXNtsNY8vHSXPxhei9/k/2XF8xOaWZWBa3ab+gEtU7Y/hCmtxpClw3bZR0Jy/9xKSzPOgoza6l8DcrMzHLJLSiADp1h8dyso2iw4ZL3YMNtsw7DzCxTTlAA2x2YdQQrWDbnmqxDWMEld83IOgQza4GcoAyAW2bCra+uuG1cehns5fcWrVS+W8d2dO/Ufi1EZmbri9pFn/LB4s/KLu8EZQCcsEPyKNRrKrw4D247dWg2QZlZizDsxtLbfZOEmZnlkhOUmZnlkhOUmZnlkhOUmZnlUlkJStIhkmZKmiVpTIn9P5L0kqQXJP1D0jaVD9WysFmHZVmHYGbrueWffjS/1PZmE5Sk1sBVwKFAH+A4SX2Kij0HDImI/sB44BdrFq7lxWYdPNaRmVXX8k8Wr16CAnYDZkXEGxHxGTAOOKqwQEQ8FBEfpatPAj3XJFgzM7NyEtRWwNsF6zXptsacAtxTaoek0yRNlTS1tra2/CjNzKzFKeeLuiqxreQMDJK+CQwB9i21PyKuBq4GGDJkiGdxaMQG7Vqz4KPyv21dLXXLlsNyd/GZWTbKSVA1wNYF6z2Bd4sLSToAOB/YNyI+rUx4LVPfLTeGjt2yDgNmtWXRkqVZR2FmLVQ5XXxTgO0k9ZbUDjgWmFhYQNIg4I/AiIh4v/JhmplZS9NsgoqIOuAM4F7gZeD2iJgh6RJJI9JiVwAdgf+TNE3SxEaqMzMzK0tZg8VGxCRgUtG2CwuWD6hwXGZm1sJ5JAkzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8ulsqbbsLWsQ2dYPDfrKGDZUkKts47CzFooJ6g82u7ArCNITL2B+Kw26yjMrIVyF5+ZmeWSE5SZmeWSE5SZmeWSr0FZk1q1ErWLPsk6DDNbny1fvqzUZicoa1Kn9m0YMXCrrMMws/XY8k8Xf1hqu7v4zMwsl8pKUJIOkTRT0ixJY0rsby/pz+n+pyT1qnSgZmbWsjSboKkvuL4AAA89SURBVCS1Bq4CDgX6AMdJ6lNU7BRgfkT8B/Ar4OeVDtTMzFqWclpQuwGzIuKNiPgMGAccVVTmKODGdHk8MFySKhemmZm1NOUkqK2AtwvWa9JtJctERB2wEOhaXJGk0yRNlTS1ttYjFOTe5v1g022zjsLMWqhy7uIr1RKK1ShDRFwNXA0wZMiQlfZbzhx6edYRmFkLVk4LqgbYumC9J/BuY2UktQE6A/+uRIBmZtYylZOgpgDbSeotqR1wLDCxqMxE4KR0+RjgwYhwC8nMzFZbs118EVEn6QzgXqA1cH1EzJB0CTA1IiYC1wE3S5pF0nI6tppBm5nZ+q+skSQiYhIwqWjbhQXLnwBfq2xoZmbWknkkCTMzyyUnKDMzyyVldS+DpEXAzEwOvnq6AR9kHcQqcLzV5Xiry/FWV97i3SYiuhdvzHI085kRMSTD468SSVMdb/U43upyvNXleKvDXXxmZpZLTlBmZpZLWSaoqzM89upwvNXleKvL8VaX462CzG6SMDMza4q7+MzMLJecoMzMLJcySVDNTSGfNUlbS3pI0suSZkj6frp9U0n3S3ot/blJ1rHWk9Ra0nOS7krXe0t6Ko31z+lAv7khqYuk8ZJeSc/zHjk/vz9MXwvTJd0mqUOezrGk6yW9L2l6wbaS51OJ36TvvxckDc5JvFekr4cXJN0hqUvBvnPTeGdKOjgP8Rbs+4mkkNQtXc/l+U23n5mewxmSflGwPdPz26iIWKsPkgFnXwe2BdoBzwN91nYczcS4BTA4Xe4EvEoy3f0vgDHp9jHAz7OOtSDmHwG3Anel67cDx6bLfwBOzzrGonhvBL6dLrcDuuT1/JJMyDkb2KDg3I7M0zkG9gEGA9MLtpU8n8BhwD0k87gNBZ7KSbwHAW3S5Z8XxNsn/T/RHuid/v9onXW86fatSQbSfhPolvPzux/wANA+Xe+Rl/Pb2COLFlQ5U8hnKiLei4hn0+VFwMsk/6QKp7a/EfhyNhGuSFJP4HDg2nRdwP7A+LRIbmIFkLQxyRvoOoCI+CwiFpDT85tqA2yQzne2IfAeOTrHETGZledga+x8HgXcFIkngS6Stlg7kSZKxRsR90UyIzfAkyRzz0ES77iI+DQiZgOzSP6PrDWNnF+AXwFns+IErbk8v8DpwOUR8Wla5v10e+bntzFZJKhyppDPDUm9gEHAU8BmEfEeJEkM6JFdZCv4X5I3yfJ0vSuwoODNnrdzvC1QC9yQdkteK2kjcnp+I+Id4ErgLZLEtBB4hnyfY2j8fK4L78FRJK0QyGm8kkYA70TE80W7chkvsD3wpbRb+hFJu6bb8xpvJgmqrOnh80BSR+AvwA8i4sOs4ylF0hHA+xHxTOHmEkXzdI7bkHQ//D4iBgFLSLqgcim9dnMUSffHlsBGwKEliubpHDcl168PSecDdcAt9ZtKFMs0XkkbAucDF5baXWJbHs5vG2ATkm7H0cDtaW9LXuPNJEGVM4V85iS1JUlOt0TEX9PNc+ub6unP9xt7/lq0FzBC0hyS7tL9SVpUXdLuKMjfOa4BaiLiqXR9PEnCyuP5BTgAmB0RtRGxFPgrsCf5PsfQ+PnM7XtQ0knAEcAJkV4gIZ/xfpHkA8vz6XuvJ/CspM3JZ7yQxPXXtOvxaZIel27kN95MElQ5U8hnKv1UcR3wckT8smBX4dT2JwF3ru3YikXEuRHRMyJ6kZzLByPiBOAh4Ji0WC5irRcR/wLelrRDumk48BI5PL+pt4ChkjZMXxv18eb2HKcaO58TgRPTu82GAgvruwKzJOkQ4BxgRER8VLBrInCspPaSegPbAU9nEWO9iHgxInpERK/0vVdDcmPVv8jp+QUmkHyARdL2JDcnfUAOz2+DLO7MILnL5VWSu0XOz/IukUbi25ukifsCMC19HEZybecfwGvpz02zjrUo7mF8fhfftiQvslnA/5HeuZOXBzAQmJqe4wkkXQ+5Pb/AxcArwHTgZpI7nnJzjoHbSK6PLSX5Z3lKY+eTpEvnqvT99yIwJCfxziK5FlL/nvtDQfnz03hnAofmId6i/XP4/C6+vJ7fdsCf0tfws8D+eTm/jT081JGZmeWSR5IwM7NccoIyM7NccoIyM7NccoIyM7NccoIyM7NccoIyM7NccoKy9ZKS6Tz+s2B9S0njm3rOKtR9kaSfpMuXSDqgQvVuoXS6lGqQtHgVyj6gHE13Yi2TE5Str7oADQkqIt6NiGOaKL9aIuLCiHigQtX9CLimQnWtqZspOH9mWXCCsvXV5cAXJU1LJ8LrVT95m6SRkiZI+puk2ZLOkPSjdGT1JyVtmpb7oqS/S3pG0qOSdiw+iKSxko5Jl+dIuljSs5JerC8vaaN0Arkp6TEam17mq8Df0+dMktQ/XX5O0oXp8qWSvp0uj07rfEHSxQUxfVPS0+nv/kdJrYti7ibpCUmHp622yWnZ6ZK+lBabCBy3muferCKcoGx9NQZ4PSIGRsToEvt3Bo4nmffmMuCjSEZWfwI4MS1zNXBmROwC/AT4XRnH/SAiBgO/T58DyTAyD0bEriSTxl2RTi/SIB0DbX6kc/UAk0mmRtiYZGTvvdLtewOPSjqIZMy03UiGjdpF0j6SdgK+AewVEQOBZcAJBcfZDLgbuDAi7k7Pwb1p2QEkQwwREfOB9pK6lvE7m1VFm+aLmK2XHopkMspFkhYCf0u3vwj0T6da2RP4v2R8WCAZf6859SPfPwN8JV0+iGTE+fqE1QH4AslEmPW2IJkjq96jwFkkM/neDRyYTvHQKyJmSjo1rfe5tHxHkoTVH9gFmJLGvQGfj2LelmRMvu9FxCPptinA9eno/RMiYlpBDO+TTC8yr4zf26zinKCspfq0YHl5wfpykvdFK5IJCQeuZr3L+Pz9JeCrETGzied9TJK46k0BhgBvAPeTTItwKkniq6/z/0XEHwsrkXQmcGNEnFviGHXp8w8GHoFk5lVJ+5DMyHyzpCsi4qa0fIc0LrNMuIvP1leLgE6r++RIJqicLelrkEzBImnAalZ3L3BmOlUHkgaVKPMq0Kvg+J+RjOz9dZLpzx8l6TJ8tKDOUWlLD0lbSepB0kI6Jl1G0qaStqmvlmSm2h0ljUn3b0My4eU1JFPMDK7/fYHNSUbpNsuEE5StlyJiHvBYeuH/itWs5gTgFEnPAzNIZtVdHZeSdK+9kN6ocWmJeJcAr0v6j4LNjwJzI5kb6VGSieQeTcvfB9wKPCHpRZJJHztFxEvABcB9kl4gaX1tUXCcZSTzhu2X3oY/DJgm6TmSmzR+nRbdBXgyPp/S3myt83QbZjkh6Whgl4i4IAex/BqYGBH/yDoWa7nWuAUl6WhJUXgLbnpL7/GrUMfjzewv+wuGZuuqiLiD/HSpTXdysqytcQtK0u0kXQj/iIiL0m3DgJ9ExBHNPLd12uXQ3DEWR0THNQrUzMzWKWvUgkov0O5FMp3wsQW7Lif5Dsc0ST8ses4wSQ9JupXklt6GFlITXxqsf27DFwzXJG4zM8u/Nb3N/MvA3yPiVUn/ljQ4Ip4l+ZJkUy2o3YCdI2J20fb6Lw1eln77fcP6HekXDCcCF0TE/WsYt5mZ5dyaXoM6DhiXLo+j/KFRni6RnCD57sfJki4C+qVfpITPv2B4tpOTmVnLsNoJKh0CZX/gWklzgNHAN+q/69GMJaU2RsRkYB/gHZIvDdYPOVP4BUMzM2sB1qQFdQxwU0RsExG9ImJrkmFZ9mY1vyTZ2JcGKfEFQzMzW7+tSYI6DrijaNtfSK4jvQDUSXq++CaJZgyj9JcGS33B0MzM1mP+oq6ZmeWShzoyM7NccoIyM7NccoIyM7NccoIyM7NccoIyM7NccoIyM7NcqliCknS9pPfTCdnqt31N0gxJyyUNqdSxqk3SD9O4p0u6TVKH5p+VLUmtJT0n6a6sY2lOqddKnjne6nK81bWuxVuoki2oscAhRdumA18BJlfwOFUlaSvgLGBIROwMtGbFkdrz6vvAy1kHUaaxrPxaybOxON5qGovjraaxrFvxNqhYgkrH0ft30baXI2JmpY6xFrUBNpDUhmRE9XczjqdJknoChwPXZh1LOUq9VvLM8VaX462udS3eQr4GVSQi3gGuBN4C3gMWRsR92UbVrP8FzgaWZx2ImVmlOEEVkbQJcBTQG9gS2EjSN7ONqnGSjiAZYPeZrGMxM6skJ6iVHQDMjojaiFgK/BXYM+OYmrIXMCKd8mQcsL+kP2UbkpnZmnOCWtlbwFBJG6ZzWw0nxzcfRMS5EdEzInqR3MzxYETktsVnZlauSt5mfhvwBLCDpBpJp0g6WlINsAdwt6R7K3W8aomIp4DxwLPAiyTn6OpMg1rPlHqtZB1TUxxvdTne6lrX4i3k6TbMzCyX3MVnZma55ARlZma55ARlZma5VO2x+K6Q9IqkFyTdIalLpY5XTZK6SBqfxv6ypD2yjqkpkjpIelrS8+kYghdnHVNzJB0iaaakWZLGZB1PcxxvdTne6lrX4m0QERV5APsAg4HpBdsOAtqkyz8Hfl6p41XzAdwIfDtdbgd0yTqmZuIV0DFdbgs8BQzNOq4m4m0NvA5sm57f54E+WcfleB2v483Xo9pj8d0XEXXp6pNAz0odr1okbUySbK8DiIjPImJBtlE1LRKL09W26SPPt2fuBsyKiDci4jOSLxgflXFMTXG81eV4q2tdi7fB2rwGNQq4Zy0eb3VtC9QCN6TTV1wraaOsg2pOOt3GNOB94P5Ivs+VV1sBbxes16Tb8srxVpfjra51Ld4GayVBSTofqANuWRvHW0NtSLoqfx8Rg4AlQO77bCNiWUQMJGml7iZp56xjaoJKbMtzi8/xVpfjra51Ld4GVU9Qkk4CjgBOiLRDNOdqgJqCFsh4koS1Tki7Ix8m3/O/1ABbF6z3JN9Tmjje6nK81bWuxdugqglK0iHAOcCIiPiomseqlIj4F/C2pB3STcOBlzIMqVmSutffISlpA5IBb1/JNqomTQG2k9RbUjuSMQQnZhxTUxxvdTne6lrX4m3QplIVpeM9DQO6pePv/RQ4F2gP3J+Mu8qTEfHdSh2zis4Ebkn/mG8AJ2ccT3O2AG6U1JrkQ8ftEZHbqd8jok7SGcC9JHcYXR8RMzIOq1GOt7ocb3Wta/EW8lh8ZmaWSx5JwszMcskJyszMcskJyszMcskJyszMcskJyszMcskJyszMcskJyszMcun/B1mhBLOQfv4RAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi2=3.396390, p=0.065339 for all events secure, exploiting aggregates\n",
      "Chi2=3.396386, p=0.065339 for all 161 time moments secure\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "from kmsurvival import main\n",
    "sys.argv[1:] = ['-i2', '--plot-curves']\n",
    "await main()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Complete Run: 5 logrank tests + survival tables\n",
    "To run the demo with three parties on localhost, for instance, we add `-M3` as command line option and run [kmsurival.py](kmsurvival.py) outside this notebook using a shell command. The plots are not shown this way, so instead we print the survival tables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using secure fixed-point numbers: SecFxp64:32\n",
      "Dataset: aml, with 3-party split, time 1 to 41 (stride 4) months\n",
      "2019-09-18 12:13:00,446 Logrank test on all events in the clear.\n",
      "Chi2=2.675902, p=0.101878 for all events in the clear\n",
      "2019-09-18 12:13:00,478 Start MPyC runtime v0.5.10\n",
      "2019-09-18 12:13:02,024 All 3 parties connected.\n",
      "2019-09-18 12:13:02,024 Logrank test on own events in the clear.\n",
      "Chi2=0.510517, p=0.474915 for own events in the clear\n",
      "          removed  observed  censored  entrance  at_risk\n",
      "event_at                                                \n",
      "0.0             0         0         0         4        4\n",
      "3.0             1         1         0         0        4\n",
      "5.0             1         1         0         0        3\n",
      "8.0             1         1         0         0        2\n",
      "12.0            1         1         0         0        1\n",
      "          removed  observed  censored  entrance  at_risk\n",
      "event_at                                                \n",
      "0.0             0         0         0         4        4\n",
      "2.0             1         1         0         0        4\n",
      "3.0             1         1         0         0        3\n",
      "7.0             1         1         0         0        2\n",
      "11.0            1         1         0         0        1\n",
      "2019-09-18 12:13:02,196 Logrank test on aggregated events in the clear.\n",
      "Chi2=2.685357, p=0.101275 for aggregated events in the clear\n",
      "          removed  observed  censored  entrance  at_risk\n",
      "event_at                                                \n",
      "0.0             0         0         0        11       11\n",
      "4.0             3         2         1         0       11\n",
      "8.0             4         3         1         0        8\n",
      "12.0            3         2         1         0        4\n",
      "44.0            1         0         1         0        1\n",
      "          removed  observed  censored  entrance  at_risk\n",
      "event_at                                                \n",
      "0.0             0         0         0        12       12\n",
      "4.0             6         5         1         0       12\n",
      "8.0             3         3         0         0        6\n",
      "12.0            3         3         0         0        3\n",
      "2019-09-18 12:13:02,337 Optimized secure logrank test on all individual events.\n",
      "2019-09-18 12:13:02,337 Interval 1 (time 1 to 4) # observed events = 7\n",
      "2019-09-18 12:13:02,337 Interval 2 (time 5 to 8) # observed events = 6\n",
      "2019-09-18 12:13:02,337 Interval 3 (time 9 to 12) # observed events = 5\n",
      "2019-09-18 12:13:02,337 Interval 4 (time 13 to 16) # observed events = 0\n",
      "2019-09-18 12:13:02,337 Interval 5 (time 17 to 20) # observed events = 0\n",
      "2019-09-18 12:13:02,337 Interval 6 (time 21 to 24) # observed events = 0\n",
      "2019-09-18 12:13:02,337 Interval 7 (time 25 to 28) # observed events = 0\n",
      "2019-09-18 12:13:02,337 Interval 8 (time 29 to 32) # observed events = 0\n",
      "2019-09-18 12:13:02,337 Interval 9 (time 33 to 36) # observed events = 0\n",
      "2019-09-18 12:13:02,337 Interval 10 (time 37 to 40) # observed events = 0\n",
      "2019-09-18 12:13:02,337 Interval 11 (time 41 to 41) # observed events = 0\n",
      "Chi2=2.675897, p=0.101878 for all events secure, exploiting aggregates\n",
      "2019-09-18 12:13:03,180 Secure logrank test for all 41 time moments.\n",
      "Chi2=2.675898, p=0.101878 for all 41 time moments secure\n",
      "2019-09-18 12:13:06,383 Stop MPyC runtime -- elapsed time: 0:00:05.905250\n"
     ]
    }
   ],
   "source": [
    "!python kmsurvival.py -M3 -i2 --print-tables --collapse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To try out other runs of the demo for yourself, remember to consult MPyC's help message, using the `-H` option:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: kmsurvival.py [-H] [-h] [-C ini] [-P addr] [-M m] [-I i] [-T t] [-B b]\n",
      "                     [--ssl] [-L l] [-K k] [--no-log] [--no-async]\n",
      "                     [--no-barrier] [--output-windows] [--output-file] [-f F]\n",
      "\n",
      "optional arguments:\n",
      "  -H, --HELP            show this help message for MPyC and exit\n",
      "  -h, --help            show kmsurvival.py help message (if any)\n",
      "\n",
      "MPyC configuration:\n",
      "  -C ini, --config ini  use ini file, defining all m parties\n",
      "  -P addr               use addr=host:port per party (repeat m times)\n",
      "  -M m                  use m local parties (and run all m, if i is not set)\n",
      "  -I i, --index i       set index of this local party to i, 0<=i<m\n",
      "  -T t, --threshold t   threshold t, 0<=t<m/2\n",
      "  -B b, --base-port b   use port number b+i for party i\n",
      "  --ssl                 enable SSL connections\n",
      "\n",
      "MPyC parameters:\n",
      "  -L l, --bit-length l  default bit length l for secure numbers\n",
      "  -K k, --sec-param k   security parameter k, leakage probability 2**-k\n",
      "  --no-log              disable logging messages\n",
      "  --no-async            disable asynchronous evaluation\n",
      "  --no-barrier          disable barriers\n",
      "\n",
      "MPyC misc:\n",
      "  --output-windows      screen output for parties i>0 (only on Windows)\n",
      "  --output-file         append output for parties i>0 to party{m}_{i}.log\n",
      "  -f F                  consume IPython's -f argument F\n"
     ]
    }
   ],
   "source": [
    "!python kmsurvival.py -H"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
